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Use of Half-Normal Plots in Interpreting 


Factorial Two-Level Experiments 


CuTHBERT DANIEL 
New York City 


Plotting the empirical cumulative distribution of the usual set of orthogonal con- 
trasts computed from a 2? experiment on a special grid may aid in its criticism and 
interpretation. Bad values, heteroscedasticity, dependence of variance on mean, and 
some types of defective randomization, all leave characteristic stigmata. The half- 
normal plot can be used to estimate the error standard deviation and to make judg- 
ments about the reality of the observed effects. An accompanying paper by A. 
Birnbaum gives some operating characteristics of these judgments. Examples are 
given of the use of half-normal plots in each of these ways. 


CoNTENTS 


. Summary 
. Construction of a half-normal grid and a sample plot 
. A test statistic for half-normal plots 
. Standardized half-normal plots 
. Some 2”* experiments 
. Use of half-normal plots in criticizing data: 
a. One defective value 
b. Two or more defective values 
c. Inadvertent plot-splitting 
d. Antilognormal distribution of error 
. Failures of half-normal plotting 
. Conclusions 
. Acknowledgements 


1. SuMMARY 


Inspection and criticism of the data from 2” experiments and from their 
fractions 2”* are always necessary. It is good to know that if the data are col- 
lected under conditions which satisfy all the assumptions of the usual factorial 
representation, then the analysis is straightforward, optimal tests are known, 
any desired significance level can be maintained, and the sensitivity of the tests 
can be calculated. But if there are one or more wild values in the set of observa- 
tions, if the error variance is not the same for all observations, or if the planned 
randomization was inadvertently modified—for example by plot splitting—then 
some means for detecting these defects should be available. 

Study of the empirical cumulative distribution of the usual contrasts com- 
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puted from the unrevised data may be useful in criticizing and interpreting 2” 
experiments. The distribution shows characteristic changes in response to each 
of the types of defect just mentioned. 

A rule of inference is given for detecting a small number of effects or inter- 
actions in the presence of a majority of “error contrasts”. This rule uses the 
ratio of the largest order-statistic of the contrasts to one whose expected value 
is nearest the standard error of the contrasts. The use of some operating char- 
acteristics (“sensitivity functions”) for this rule, developed by A. Birnbaum 
(2), is illustrated by numerical values for (p — g) = 4, 5, 6, 7 and for the case 
of only one real effect. 


2. CoNSTRUCTION OF A HALF-NORMAL GRID AND A SAMPLE PLOT 


The 31 contrast-sums from a 2°° experiment on penicillin production, de- 
scribed in Davies (7, Example 9.2), are arranged in decreasing order of absolute 
magnitude in Table 1. The largest values indicate the effects most likely to be 
non-zero. We need, of course, some rule of inference that will help us to be 
objective in judging whether or not the largest effects are real. 

This becomes a rather important matter in large factorials since in these the 
largest contrast, even in null experiments, has an expected value of several 
times the average value of the set. Thus in a 32 run experiment the expected 
value of this ratio is roughly 2.4. Using the usual ‘0.05 level of significance” in 
such experiments will result in some effect being wrongly judged real in over 
half of all experiments done. If this fact is recognized and accepted, no harm 
is done. A ready way to keep the frequency of false positives per experiment at 
a desired level will be given in the next section. 


TABLE 1 


Ordered Contrast-sums from a 2° Factorial Experiment 
(See reference 7) 


Jalues multiplied by 100 and rounded 
Value 


E 

A 

Cc 

CE 

ABCDE (Blocks) 
AB 

ABCD 
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If there were no effects, then the ordered absolute values would have the 
distribution of the range in samples of size two. The distribution of the range of 
pairs from the normal distribution is known and tabulated. Even if the distri- 
bution of error is not normal, we might expect the distribution of contrasts 
(which are sums) to approach normality quite closely. However, some inter- 
esting exceptions are given in section 6 below. 

If there is one large effect, we can expect the observed distribution of the 
ranked contrasts to be only slightly disturbed since only one error-manifesta- 
tion (that attached to the large effect) has been abstracted from the set. We 
might then expect the remaining (n — 1) contrasts to behave like a sample 
from the null population. 

If, finally, a considerable proportion of the n effects differs appreciably from 
zero, then (in the absence of knowledge of the magnitude of cz , the run-to-run 
or plot standard deviation) little can be learned from inspection of the ordered 
contrasts. In such cases, the safeguard of an independent estimate of error is 
obviously mandatory. The question of what proportion of large contrasts is a 
considerable one is partly answered in this section. 

My former practice of comparing the sample empirical cumulative distribu- 
tion (ecd) with the standard tabled distribution by plotting both on normal 
probability paper, was simplified by A. Birnbaum. He pointed out that this is 
the “half-normal distribution” with density: 


f(a) = — ere’ = for x > 0 
0, for x < 0. 
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‘Fig. 1—General purpose half-normal grid. 
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As first suggested in (6), a “half-normal probability paper” (on which the 
graph of the theoretical cumulative distribution is a straight line through the 
origin) can easily be prepared. On a sheet of “arithmetic (normal) probability 
paper” delete the printed probability scale P, for the range of P less than 50%. 
For the range of P greater than 50% replace each value of P by the corre- 
sponding value P’ = 2P — 100. A rough ruling for general purposes is given in 
Figure 1. The relation 


P’ = (i — 9)/n; += 1,2,3,+--n, 


is used for plotting the empirical distribution. The abscissae are the absolute 
values of the contrasts. 

Special-purpose grids for n = 15, 31, 63 and 127 are simpler to use since no 
probabilities need be computed. The rank numbers are given directly on the 
grid as shown in Figure 2 a, b, c, d. With these it is only necessary to piot the 
absolute values of the ordered contrasts directly on the horizontal lines of the 
grid. These grids are precise enough for all the uses proposed in this paper. 

This preparation and use of “half-normal probability paper” is illustrated in 
Figure 3, with the experimental data of Table 1. The largest contrast, that for 
E, is plotted on the line numbered 31, and so on down to the next-to-smallest 
(magnitude 0.025, name AE) which is plotted on the line numbered 2. There is 
little advantage in plotting each contrast of the smaller half of the set because 
of the enforced close correlations in magnitude between adjacent ordered values. 

Even in advance of detailed information about the sampling fluctuations of 
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Fic. 2—Half-normal scales for (a) 15 d.f., (b) 31 d.f., (c) 63 d.f., (d) 127 d.f. 





USE OF HALF-NORMAL PLOTS 


| E a 


HAA HEE 


| | Ha A CE -CE| | SERRGREE 
LT TT TR APEC 





° 


Fia. See plot of a 2° experiment. Dots and solid line-31 contrasts. Crosses and 
dashed line-24 smallest contrasts. 


order statistics, we would want to judge A, C, and E£ significantly large. Both 
on the grounds of their simple names, as well as because of their magnitudes 
relative to the remainder of the set, these three should be judged real. 

If the experiment had been a null one—all factors having no effects—the 
standard error of the contrast-sums, ¢, , could be roughly estimated by that 
contrast for which P’ is most nearly 0.683, or 68.3%. This is the 22nd largest 
for n = 31, since (22 — 4)/31 = 0.693. But this is clearly not a null experiment 
and so we ‘will be using fewer than 31 contrasts to estimate error. The value 
0.53 should not be taken to estimate ¢, . 

We are told (on page 418 of reference 7) that CE was judged likely to be 
appreciable from previous information. Similarly, week-to-week variations were 
known to be common and so this experiment was done ‘in two blocks, con- 
founding the interaction ABCDE with weeks. Since just these two contrasts 
appear next in order of magnitude, it seems fair to remove them from the set 
to be used to estimate error. 

The next two contrasts measure AB and ABCD. These and all the smaller 
ones will be (still tentatively) taken as error contrasts. We can either plot the 
26 on a grid like Figure 1, or we can move each of the 26 upward on Figure 3 
to a point corresponding to the new P’ = (i — 4)/26. It is noticeable that even 
after the removal of five contrasts, the largest one remaining is still plotted 
nearly on the old topmost line (at 0.981 versus the former 0.984). The revised 
plot is indicated on Figure 3 by the points x. The standard error of the con- 
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Fic. 4—Half-normal plots of ten sets of 31 random standard normal deviates. 
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trast sums is now estimated as 0.33. The upper end of the plot suggests, but 
hardly proves, that some truncation has occurred. Nothing will be greatly 
changed by deciding to include or exclude the block contrast from the error 
system. We can, then, safely use the whole set of 27 to estimate error. 

An estimate, ¢, , of the run error can be obtained from the obvious formula: 


ép = 6./Vn+1. 


The estimate just made, based on 27 degrees of freedom, agrees almost exactly 
with that given in reference 7, estimated there from the 15 degrees of freedom 
associated with the three-and four-factor interactions. 

A. step toward greater objectivity in deciding what is random and what is 
systematic in unreplicated factorial experiments has been made by J. W. Tukey 
(16), and has also been put forward by S. C. Pearce (12). This is to “nominate” 
all effects, interactions and block contrasts that are thought likely to be im- 
portant, before the experiment is completed. The corresponding contrasts are 
then to be excluded from further arbitrary allocation to error or to “reality’’. 
Only those not nominated will be studied for possible allocation to error. This 
suggestion seems to me a most acceptable compromise between the risky 
practice of calling everything error that looks like error, and the conservative 
practice of calling nothing error that is not found by strict randomization of 
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Fic. 5—J. W. Tukey’s half-normal grid. “For (logarithmically) plotting magnitudes of the 
2" — 1 contrasts in a 2” or 2"/2"-™ experiment. Plot every 2"~*-rd value on solid verticals; 
add tail values; add intermediate values to taste.’’ Dots are from Yates’ 25 experiment on 
beans. (17). 
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identical treatments. It is worth pointing out, however, that this rule would 
not detectibly change any conclusions or estimates made in the experiment 
just discussed. 

To give some idea of the sampling variation inherent in plots like that of 
Figure 3, Figure 4 shows half-normal plots for ten sets of 31 random standard 
normal deviates (rsnd) taken from Dixon and Massey’s tables (8). The first 31 
in each of the columns 11 to 20 are used. The printed abscissa values are correct 
only for the first plot in each set. Each sloping line should start at zero. Each 
has the slope required by the population standard deviation (1.00).and not the 
slope of the best line through the points shown. 

J. W. Tukey (16) prefers a grid that uses the logarithms of the absolute values 
of the contrasts for one coordinate. It is his judgment that the relative magni- 
tudes of the contrasts are intuitively important, and that the corresponding 
expansion of the scale for the smaller contrasts is useful in looking for some 
types of pathological error distribution. He also inclines the contrast-axis at 
135° to the horizontal positive direction instead of at right angles as proposed 
above. The purpose of this inclination is to bring the expected position of the 
plot for a null experiment to a horizontal straight line, instead of having it 
appear as an inclined line as in the grids shown in Figures 3 and 4. Figure 5 
shows this grid, sketched by the writer with permission, with Yates’ 2° experi- 
ment on beans (17) plotted as an example. The directions on the grid 
are Tukey’s. 


3. A Test Sratistic For HauF-NorMAL PLots | 
The usefulness of half-normal plots in distinguishing between the cases of no 


effects and of one effect is of special importance. This has been investigated by 
A. Birnbaum (2). The distributions of interest are those of 


Un 
Ua 


(A) 
eS = 


where: 
u, is the largest in absolute magnitude of n observed values of a random 
normal variable with population mean zero, variance o”, but with one value 
chosen at random changed by the addition to its signed value of Ac; 


u, is the absolute value of the order-statistic from the same set of » that is 
numbered nearest (0.683n + 0.5), and which therefore most nearly estimates 
o directly,in a null experiment. For values of n of 15, 31, 63, 127, the cerre- 
sponding o-estimators are the 11th, 22nd, 44th and 88th order statistics; 


i‘ is then the ratio of the largest in absolute value of the n observed values 
(one perturbed by the addition of Ac) to the corresponding o-estimator, 4, . 


A. Birnbaum has studied the distributions of ¢S“’ for n = 31, 63, 127, and 
for A’s of 0 to 6. 

A machine sampling of 2500 sets of 31 rsnd was carried out at his request 
by G. Lieberman and associates at Stanford University. This permitted an esti- 
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Fig. 6—Empirical cumulative distribution of logio ts:“, plotted on normal probability paper. 


mation of the distribution of ¢{! with good accuracy over the range of interest. 
The distribution of logyts}’ is quite closely approximated by a normal distri- 
bution with mean 0.35 and standard deviation 0.11. Figure 6 gives at the circled 
points several values from the empirical distribution of logarithms plotted on a 
normal probability scale. H. M. Truax (15) making slightly different assump- 
tions found values closely approximating those given above. 

Birnbaum’s paper gives (asymptotic) approximations for the distributions of 
i? and uses them to approximate the distributions for n = 63 and 127. The 
writer has estimated the distributions of ¢{f? using 198 samples of 15 rsnd taken 
from Dixon and Massey’s tables (8). The cumulative distributions of logiot{” 
for n of 15, 31, 63 and 127 are all shown on a normal probability scale in Figure 
7, and with comparable numerical accuracy in Table A below. 


Taste A 
Approximate Distribution of logio( ta). 
All are nearly normal for P > 0.1 with estimated parameters A, ¢. 
n a é 
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31 
63 








CUTHBERT DANIEL 


ia ae 

| PPE 
ee ee 
aa ede id 


80 85 
Per cent 


Fic. 7—Estimated cumulative distributions of logio é,“ for n = 15, 31, 63, 127. 


The relative frequency with which ¢, (the ratio u,/u, obtained as a statistic 
from an experiment) will detect a single effect of magnitude A in sets of inde- 
pendent rsnd (when the frequency of false detection per null set is set at a) is 
called the sensitivity function and is symbolized by y(n, a, A). Upper and lower 
bounds on this function found (2) for n = 31, a = 0.05, 0.2 and 0.4 are shown 
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Fic. 8—Bounds on sensitivity function of t3:“ for one effect of size A (after Birnbaum). 
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Fie. 9—Approximate sensitivity function of tis‘ for one effect (based on 198 sets of random 
normal deviates). 


graphically in Figure 8, by the three pairs of lines converging to the right. 
Empirical sampling by the writer for A = 4, using 99 sets of 31 rsnd gave the 
values shown by the three solid dots in that figure. 

A similar graph for n = 15, estimated by the writer from a sample of 198 sets, 
is given as Figure 9. Since A’s of 2 and 4 only were used, the straight lines drawn 
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are entirely conjectural, based on the analogy with the linear beunds found by 
Birnbaum for n = 31. 

The distributions of ¢,°) , °° , ete., are also of interest since more than one 
non-zero effect is usually hoped for. This has been examined empirically by 
the writer with the same 198 sets of 15 discussed above. The resulting empirical 
cumulative distributions are shown in Figure 10. 

Since the values of ¢S*’ and-of y(n, a, A) were computed for rsnd, a scale 
correction must be made when Figures 8 and 9 are used to estimate the actual 
sensitivities of factorial experiments in detecting average effects. The A of 
Figure 8 must be divided by 24/2, and that of Figure 9 by 2 to get estimates 
of sensitivity in terms of the run-to-run standard deviation. 


4. STANDARDIZED HatF-NorMAL PLotTs 


It may be convenient to use a half-normal grid with fixed guardrails for 
studying the results of a 2* experiment. Dividing the computed ranked con- 
trasts by u,, will give a scale-free set of order-statistics that should in the absence 
of real effects fall around the standard line drawn in Figure 11a. From the ecd 
of logy ¢{3? given in Table A and in Figure 7, it is easy to compute values of 
t{? that will only be exceeded in, say, one experiment in twenty. Thus at prob- 
ability 0.05 we read from Figure 7 a logarithm of 0.49, corresponding to a ¢,; of 
3.09. The expected value of this ¢ is 2.13, and so the excess of 0.96 gives the 
distance of the 0.95 guardrail from the standardized expected position. Table 2 
gives these critical deviations for the usual four values of n, and for error-rates 
of 0.05, 0.2 and 0.4 false positives per experiment. 

The three values in the first line of Table 2 places the offsets of the three 


dotted lines in Figure lla at ordinate 15. The corresponding critical deviations 
for the next two smaller standardized order-statistics, deduced from Figure 10, 
lie approximately on straight lines as indicated in Figure 11a. 

In factor-screening experiments it seems desirable to use values of a of 0.4 
or even 0.8, knowing that this will produce a number of “false alarms”. The 


TABLE 2 
Critical Deviations for Standardized Half-normal Plots 


15 
31 
63 
127 


nm = number of orthogonal contrasts 

r = expected value of ¢, 

a = frequency of false positives per experiment 

d = cell entries, deviations to right of standard line. 
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justification is of course in the increased sensitivity to small effects as indicated 
in Figures 8 and 9. Erroneously selected effects have good prospects of being 
exposed in later work. Missed real effects on the other hand are likely to be 
dropped from study, and if they are numerous, loss of knowledge of their com- 
bined effects may be serious. 

Standardized half-normal grids for n = 31, 63 and 127 are shown in Figure 
11, panels b, c, d. 


5. Some 2?-* EXPERIMENTS. 


Figure 12 shows the appearance of eight actual experiments, with p-q equal 
to 5, on a standardized half-normal grid. 


. Yates’ 2°-° on mangolds (reference 11, page 269) 

. Bennett and Franklin’s 2°™* (reference 1, page 586) 

. Davies’ 2°~° on penicillin (reference 7, page 417) 

. Cochran and Cox’s 2°~* (reference 4, page 254) 

. Daniel and Riblett’s 2°~* (reference 5) 

. U.S. Steel Corp., Applied Research Laboratories oe 

. General Foods Corp., Research Center 2**~° 

. U.S. Steel Corp., Applied Research Laboratories 2°°-*°. 


To save space four plots have been made on each grid. The indicated abscissa 


TaBiz 3 
Standardized Ordered Contrasts for Eight 2° Experiments 


(See Figure 12) 
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Fic. 12—Standardized half-normal plots of eight 2?-¢ experiments; p — q = 5. 
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TABLE 4 
Standardized Ordered Contrasts for Six 2* Experiments 


(See Figure 13) 


c 
24-0 


1.19 
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01 
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18 
.23 
45 
51 
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* not un but roughly us because of three large effects 


values should of course be decreased by / for each successive plot. The standard- 
ized contrasts are given in Table 3. 

Each of these experiments merits more discussion. For example, there appears 
to be a defective value (specified in the original paper as d) in the experiment 
graphed in plot a. Revision of this value smooths the graph, gives an estimate 
of the error standard deviation with 27 df closely matching that given by 
Kempthorne (11, p. 270) with 13 df, and increases slightly the apparent signifi- 
cance of the four largest effects. These discussions might take other directions 
as indicated later in section 6. 

Figure 13 gives six examples of 2* experiments on the appropriate grid. The 
corresponding standardized ordered contrasts are shown in Table 4. 


. U.S. Steel Corp., A. R. L. 2*~° 

. Okonite Corp. 2*~° 

. Davies, O. L. 2*~° (reference 7, page 276) 

. Davies, O. L. 2°~* (ibid. Page 462) 

. General Foods Corp., Research Center. 2°~* 
. Ibid. 2”~* 


Figure 14 shows three 64-run experiments 


a. General Foods Corp., Research Center. 4” X 2° in half replication. 

b. U. S. Steel Corp., Applied Research Laboratories. 47 X 2* in quarter 
replication. 

c. Ibid. 4” X 2° in eighth replication. 
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Fie. 14—Standardized half-normal plots of three 2?-¢ experiments; p — q = 6. 
Fie. 15—Standardized half-normal plots of two 2?-¢ experiments; p — q = 7. 
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Fia. 16—Effect on half-normal plot of one defective observation. (a) contrasts from original 
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Figure 15 gives two 128 run experiments. 


a. General Foods Research Center. 4’ X 2° in eighth replication. 
b. Ibid. 4 X 2° in sixteenth replication. 


The experiments were chosen as fairly ‘“good”’ examples without the compli- 
cations to be mentioned below. 


6. Use or Hatr-NorMAt Puiots in Criticizine Data. 


The half-normal plot appears to be sensitive to certain types of defects in 
experimental data and can sometimes be useful in detecting these defects. 
Some experience has been gained with four types of defects: 


a. One defective value 

b. Two or more defective values 

c. Inadvertent plot-splitting 

d. Antilognormal distribution of error. 


a. One defective value 


A single very wild value—too large or too small by, say 5 og—is easily spotted 
by inspection of the raw data if only one or two factors have noticeable effects. 
But if more factors (or interactions) are influential, and if the wild value is only 
off by 2 og , then there may be several responses exceeding the maverick and 
thus masking it. Furthermore most of the contrasts used for error will be in- 
flated by the maverick so that the real effects may not be detected. 

The mistaken value will appear on one side of every contrast and will on the 
average increase its absolute value. Panel a of Figure 16 shows the half-normal 
plot for the data given in Table 5. The experiment was a 4” X 2? in quarter 
replication. Factors B and C and their interaction were used. to represent a 
four-level: quantitative factor. Factors D and F, with DF, were assigned to a 
four-level qualitative variable. ‘The alias sub-group contained ABDE, ACDF, 
BCEF. © 

The visually obvious failure of the smaller contrasts to “point toward the 
origin suggests an error of about 7 in all the contrasts. Inspection of the data 
makes one think that measurement beef is at fault. Revising this value from 
28.7 to 21.7 changes each contrast by + 7. It is not necessary to recalculate the 
contrasts from the data. Table M of reference (7) shows the sign of the revision 
required for each treatment combination and for every effect, 

A half-normal plot of the revised contrasts is given in Figure 16 panel b. Now 
selecting B, BC, and C as real effects, the remaining 12 contrasts give the error- 
plot shown by the x-points and by the dashed line in the same panel. One might 
guess that’ one or both of the contrasts AB and BD also measured real effects 
but the data do not suffice to decide this matter. 

Bad values, whose magnitudes are due neither to the regular extor-eystem. 
nor to the effects of the deliberately varied factors, are a major hazard in 2”* 
factorials. ‘They are more damaging to the experiment than missing values 
because: they are not easily identifiable. Nor are such values as rare as one might 
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TaBe 5 
Data and Contrasts from a Defective (4° X 2*)/4 
(See Figure 16) 
B, C, BC used for a quantitative, evenly spaced four-level factor 
D, F, DF used for a qualitative four-level factor 


Experimental Unadjusted Unadjusted Adjusted* 
Conditions Data Contrasts Contrasts 


(1) 
aef 
be 
abf 
cf 


ace 
beef* 
abc 
def 
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bdf 
abde 
cde 
acdf 
bed 
abcdef 


290.5 
—.5 
42.9 
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* beef was adjusted by —7.0 
Alias subgroup I, ABDE, ACDF, BCEF. 


hope. If my own experience is typical, they appear with relative frequencies 
from 0.1 to 0.01 depending on the complexity of the experimental situation 
and on the experience of the experimenter. 

Many research workers feel that 2 f. 7. c. 2””* experiments must fail in their 
work because of the known presence of higher order interactions, i.e. of re- 
sponse surfaces highly curved in the region of interest. In my opinion a large 
part of this feeling is based on the appearance of erroneous values which were 
not identified as such. 

A common feature of experiments with an evenly spaced four-level quantita- 
tive factor is the appearance of both the “linear” contrast (here B) and of the 
“cubic” contrast (here BC). The latter often appears with reversed sign at 
about half the magnitude of the former. The “linear” contrast actually measures 
the linear effect minus half the cubic component, while the “cubic” contrast 
measures the cubic effect plus half the linear. A large linear effect will then 
influence both contrasts. The most reasonable interpretation for the current 
experiment is that the same linear component is producing a B-contrast of 43 
and a BC contrast of — 20. 


b. Two or more defective values 


Two mavericks with roughly the same magnitude of bias will act together 
to inflate the absolute values of half the contrasts in a 2” experiment. The 
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shape of Figure 13f suggests that some such thing has happened. If the two do 
not have the same bias then all the contrasts will be disturbed, half of them by 
the difference in bias of the two bad values, half by the sum of the two biases. 


Similarly if three mavericks occur, one-fourth of the contrasts receive the 
maximum disturbance. 


c. Inadvertent plot-splitting 


Table 6 gives the contrasts obtained in a 2°~° designed to study the effects 
of five composition variables and of one process variable on the tensile strength 
of a certain type of steel. 

Metallurgical common sense dictated the splitting of each of the 32 heats 
(each of a different chemical composition) into two parts, to accommodate the 


TABLE 6 
Standardized Ordered Contrasts from a Split-Plot 2°-° 


(See Figure 17) 
Within Plot Whole Plot 


13 
23 
42 
49 
58 
77 
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Fia. 17—Half-normal plots of contrasts from a 2° split plot experiment. (a) combined, 
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two processing conditions to be used. This natural way of doing the experiment 
corresponds of course to a split-plot design. If all four and five factor interaction 
contrasts are used for error, a large number of 2 and 3fi appear “highly signifi- 
cant”, due entirely to the predominance of within-plot contrasts in the terms 
used for error. 

The half-normal plot for the whole collection of 63 contrasts is shown in 
Figure 17a. The shape of this plot is due to the mixture of two error variances. 
Splitting the contrasts into those containing F (the process variable and hence 
the plot-splitter) and those containing only A, B, C, D and E, gives the two 
plots shown in Figure 17b and c. (For convenience, since F was clearly influ- 
ential, the F contrast was omitted in making this figure. An ordinary 31-con- 
trast grid could then be used.) The within-plot standard deviation appears to 
be about one-fourth that for whole-plots. Only three 2fi and no 3fi are now 
judged real.. 

The half-normal graph was used here not so much to criticize the data as to 
criticize the’ statistician’s defective understanding of the experiment, especially 
of its randomization scheme. 

The appearance of Figure 17c again raises the suspicion that two maverick 
heats were present. Half the contrasts,appear to be “too large” for the pattern 
set by the 16 smallest ones. This and other problems a with this ex- 
periment will be‘dealt with elsewhere. (10) 


d. Antilogtiormal distribution of error 


When there are reasons for thinking that the error of runs increases with the 
population mean response, it is common practice to try to standardize the data 
by using the logarithms of the observations for analysis. A useful discussion of 
the limitations of this transformation is given by Severo and Olds (14). A distri- 
bution of random numbers whose logarithms are normally distributed is called 
logarithmico-normal by them; the term antilognormal is used below. 

A sample from an antilognormal distribution can be synthesized by taking 
the antilogs of a set of random normal numbers. Putting y. = In y, , where 
Yo is N(us , 03), and 4, is antilognormal (; , 07), we have that m, the coefficient 
of variation of y, is related to o3 by 


= In(m’ + 1). 
If o} is set equal to 1, then 


m= Ve" 1 = Ve—1= 131. 


A set of rsnd taken from Column 77 of reference (8) were so transformed and 
then put through the standard computation as if they were the data from a 2° 
experiment. ‘The resulting contrasts are given in column I of Table 7. The half- 
normal plot is shown in Figure 18a. The sag observed is characteristic although 
with real data the sag may be partly or entirely concealed by the presence of 
real effects. The smaller two-thirds of the contrasts in a 2° will usually still 
show the sag if the coefficient of variation is as large as 1. 

Mathematical analysis by H. Scheffé (13) shows that the expected shape of 
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TABLE 7 
Effects of Antilognormal Error on the Half-normal Plot 
(See Figure 18) 
Column I, Contrasts from a synthetic set of 32 numbers, CV = 1.31 (18a) 

II. Data from a 25-1 
III. Contrasts from unlogged data (18b) 

IV. Logged data 

V. Contrasts from logged data (18c) 


II iil 


2.2 


3.5 


the plot for antilognormal data is always convex downward, and more so for 
larger coefficients of variation. It does not follow that an observed sag always 
implies antilognormality. For example, a single maverick may produce the 
same appearance. In my opinion, the half-normal plot is not sensitive to and 
hence not very useful in detecting antilognormality. 

The example given in Table 7 (columns II-V) and in Figure 18b and c, sug- 
gests that logging may be advisable even if the plot shows no sag. Figure 18b 
gives the half-normal plot made by the writer from the data from a 2°"* experi- 
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Fic. 18—Effects of antilognormal error on half-normal plots. (a) Contrasts from synthetic 
data with coefficient of variation of 1.31; (b) Contrasts from a 25 with coefficient of variation 
0.6; (c) Contrasts from logarithms of data of b. 
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ment given in standard order in Table 7 column II. Here all 15 contrasts fit 
neatly on one line. The run error deduced from this plot is 195/+/16 or 49. 
The mean of the responses, which are intrinsically positive, is 42. The estimated 
coefficient of variation of 1.2 suggests logging. A less casual examination of the 
data would have revealed that numbers varying from 2 to 191 might well not 
all have the same precision. 

Using the logarithms of the numbers, as given in column IV, the contrasts 
of column V were computed and plotted in Figure 18c. It now appears that at 
least four and possibly five of the factors varied had consistent multiplicative 
effects. The coefficient of variation is now estimated at 0.6. 


7. FarLurREs oF HAtF-NorMAL PLOTTING 


It is unnecessary to warn experienced statisticians that the use of half- 
normal plots suggested here is still full of subjective biases, that it is not offered 
as a general substitute for the analysis of variance, and that its use in a routine 
way may be catastrophic. More optimistic and less experienced statisticians 
may get the impression from the successful examples given that a panacea is 
being offered that can hardly fail. This is not the case. 

The writer has made half-normal plots for all the 2”* experiments given in 
the standard texts and in the easily accessible journals. Over nine tenths of 
them give tolerably linear graphs for the smaller contrasts. A similar proportion 
holds in several hundred unpublished industrial experiments. No generic dif- 
ferences between published and unpublished experiments have been noticed. 

Smoothness of graph is not a guarantee of perfection, as the example just 
given in section 6d shows. Errors of both sorts may be expected to occur. A 
smooth line (and no effects securely found) might sometimes lead us to over- 
estimate the magnitude of the error and hence to miss a number of real effects. 
An extremely irregular line might lead to the other sort of error, judging effects 
to be real when in fact they are not large. Wide variation in error from run to 
run introduces correlation among the contrasts. This will often lead to an ir- 
regular graph with bunched points. No quantitative results are given here, but 
it must be reported that some few experiments have shown very irregular 
graphs. 

The experimenter who knows his run-to-run error quite exactly, who knows 
which are his influential factors and what their influential ranges, may well 
design and carry through a 2” experiment in which nearly all main effects 
and many interactions are large compared to their errors. In such situations, 
half-normal plots will be of little use. 

The 2°-* described by Box (3) was used to derive estimates of 15 constants 
in a second-order incomplete quadratic equation, the five pure quadratic terms 
being confounded with the mean in the first phase of this work. The correspond- 
ing 15 contrasts are plotted in Figure 19. The run error, known from later work, 
would lead us to expect the dashed line in that figure. The 15 effects appear to 
fall into a nearly normal distribution with their own variance, much larger than 
the error variance. The reality of many of these terms is proved by later experi- 
mental work. 
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Fie. 19—Half-normal plot of a 25! experiment with many effects. Dashed line shows true error. 


There may be some sort of central limit theorem analog operating here on 
the (discrete) population of effects, but no success can be reported in its formu- 
lation. 


8. CoNCLUSIONS 


The use of a plot of the empirical cumulative distribution of the usual orthog- 
onal contrasts from 2” * experiments (for p — g from 4 to 7) has been described. 
It has some value as an indicator of certain common types of defect in the 
data from industrial experiments. When only a small proportion of the totality 
of contrasts have effects, this plot can be used to make judgments about the 
reality of the largest effects found. 

Three sets of questions arise in the writer’s mind. None of them is answered 
here. The first set concerns obvious variants of the method of inference pro- 
posed in this paper and in the accompanying one by A. Birnbaum. Should an 
invariable rule be set up, using only some fixed proportion of the smaller con- 
trasts to estimate error? Should one use the “error contrasts” to form a mean 
square error term? Should one decline to use only higher-order interactions for 
error since some plot-splitting is almost inevitable in multi-stage processes? 
Should one insist on at least partial duplication of 2‘ experiments when no 
good previous estimate of error is available? 

The second set of questions concerns the consequences of the inveterate habit 
of studying parts of a statistical procedure, and of then using the conclusions 
found for the parts as if this guaranteed that the entire procedure must now 
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somehow be better. Like most statisticians, I am in the habit of using several 
graphical methods in the studying the results of a 2”* experiment. Two of the 
most useful ones employ the residuals from the regression equation implied by 
the significance judgments made. The ecd of the residuals is sometimes useful 
in spotting one or more extreme discrepances. A plot of the residuals against 
their respective regression values, as suggested by Anscombe and Tukey, has 
often been very revealing, especially as an indicator of some dependence of 
variance on mean value. The combined use of all three practices has not been 
studied, and yet it has been implicitly assumed that the combination has de- 
sirable operating characteristics against the alternatives most commonly met, 
and at some desirable overall significance level. 

The third set of questions plaguing the consience of one applied statistician 
concerns the repeated use of univariate statistical methods instead of a single 
multivariate system. I know of no industrial product that has but one important 
property. And yet, mainly because of convenience, partly because of ignorance, 
and perhaps partly because of lack of full development of methods applicable 
to industrial research problems, I find myself using one-at-a-time methods on 
responses, even at the same time that I derogate a one-at-a-time approach to 
the factors influencing a response. To what extent does this simplification 


invalidate the judgments I make concerning the effects of all factors on all 
responses? 
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On the Analysis of Factorial Experiments 
Without Replication’ 


ALLAN BirNBAUM 
New York University 


This paper discusses several methods for judging which of m contrasts provided 
by a factorial design without replication may be different from zero. These include 
the technique of half-normal plotting, proposed by C. Daniel, for which some operating 
characteristics are given. 


1. SuMMARY 


Some questions of validity of formal models and statistical inferences are 
discussed briefly as they arise in connection with factorial experiments without 
replication in exploratory research. It is suggested that under certain conditions 
(which are reviewed) the following schematic model is useful: The m contrasts 
(estimates of effect-parameters) y; are independently normally distributed with 
common unknown variance; at most (any) r of their means differ from zero 
(here r is a specified integer, 1 < r < m). The relevant problem is to infer, on 
the basis of one set of observed values y; , --* , Ym, Which, if any, of the means 
are nonzero. Optimal inference rules are described for r=1 (in which case they 
are shown to depend on the statistic max,y?/)., y?) and r = 2. An alternative 
graphical procedure developed by C. Daniel [1] is related to a statistic ¢,, which 
may be regarded as an approximation to the optimal statistic for the case 
r = 1. Some important advantages of procedures based on ¢,, are discussed; 
critical values, power, and related properties of such procedures are given; 
and comparisons with more familiar statistical procedures are made. 


2. INTRODUCTION 


The usual methods of analysis and interpretation of data from factorial 
experiments without replication (analysis of variance, Model I, and related 
methods) are based on certain well-known theoretical assumptions. The validity 
of these methods depends most sensitively on the assumption that certain 
high-order interactions are zero—this assumption is crucial especially in the 
case of fractional replicates. 

Their validity is moderately sensitive to the assumptions of equality of 
variances and of independence, and rather insensitive to non-normality. 


1 This research was sponsored in part by the Office of Naval Research under Contract 
Number NONR-266(33), Project Number NR 042-034, at Columbia University. Reproduction 
in whole or in part is permitted for any purpose of the United States Government. 
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One important use of these methods is made in the early exploratory stages 
of many research investigations. It is inherent in such situations that any 
statistical model, which provides a formal basis for drawing informative infer- 
ences from limited experimental data, will typically be too schematic, too 
restrictive of unknown conditions, to be claimed “valid”, or a basis for infer- 
ences which are “valid” except in the hypothetical formal sense. Such a model 
may be, perhaps along with alternative somewhat different models, a basis for 
what may be termed “plausible inferences’, that is, inferences drawn in a 
formally valid manner, based on a model which is more or less plausible (on the 
basis of previously available knowledge and/or of the experimental data under 
consideration). 

For such applications, there seems to be room for methods of analysis of 
factorial experiments without replication alternative to or complementary to 
the standard methods. This paper, and the accompanying paper [1] by Daniel, 
present respectively the statistical theory and some representative applications 
of some new methods for such experiments. The following section presents the 
statistical model which is the formal basis for these methods, and an account 
of considerations which in many experimental situations seem to make the 
model a usefully plausible one. 


3. AN INFERENCE PROBLEM RELATED TO THE ANALYSIS OF FACTORIAL 
EXPERIMENTS WITHOUT REPLICATION 


In one replication of a 2” experiment, with the usual assumptions of the 
analysis of variance, Model I, a familiar transformation of the N = 2” obser- 
vations gives N independent normally distributed statistics y, , --- , yw , having 
a common unknown variance o” (equal to the “error variance” of one observation 
divided by N), and with respective expected values equal to the “grand mean” 


up; the main effects of the p factors, a“**, ---, a**; the (?) two-factor int- 


2. 
actions, oo. mee, aAe-24. tees and the p-facto interaction qgttersedn 


The familiar standard inference methods here require the assumption that 
certain of the a’s are zero; then the sum of corresponding y7’s is used to estimate 
o, and ¢ or F statistics (or others) are used for inferences about remaining 
a’s. Similar remarks apply to fractional replicates of 2” experiments. In ex- 
periments with factors at more than two levels, an appropriate transformation 
of the observations gives N independent statistics y, , --- yw ; the main effects 
of a factor with J levels are represented by a certain (J — 1) y,’s, etc. For the 
sake of definiteness and simplicity, the following discussion of methods is given 
in terms of single replicates of 2” experiments; the simple modifications required 
for applications to other experiments without replication are left to the reader. 

Interest in an alternative to the standard assumption that certain high-order 
interactions are zero stems from the following considerations: 

1. Experience in certain areas of experimentation suggests that usually the 
number of a’s actually differing appreciably from zero in any one experiment 
will be quite small, relative to N; in particular, that it will be much smaller 
than the number of a’s which would usually not be assumed zero a priori in the 
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standard methods of inference. E.g., if all 3-factor and higher order interactions 
could be ‘assumed zero, there would remain 


r+) Ma 


statistics yo , *** , Yinc+1)/21+1 Which might be investigated by the standard 
methods individually or in groups just for inferences about the main effect and 
two-factor interaction a’s. The suggestion above implies that stronger inferences 
might be obtained by using the additional a priori assumption that, at most a 
certain small number among the statistics yz , -*+ , Ytec+1y/2}41 have non-zero 
means. 

2. Experience in certain areas of experimentation also suggests that high- 
order interactions are sometimes actually different from zero. Hence it is desirable 
to avoid the a priori assumption underlying the standard method, that all 
k-factor through n-factor interactions are zero, for some k > 1, k S p. 

The preceding points suggest as a basis for inference a new underlying assump- 
tion like the following: Any among the 2” — 1 statistics y, , --- , yy may have 
non-zero means, but at most a certain small number r of their means differ from 
zero. The statistical problem, then, is to infer which, if any, of the a’s differ 
from zero. 

In exploratory research situations, both the underlying assumption of the 
standard method and that of the preceding paragraph seem somewhat schematic. 
The method presented in the following Section 4 has a graphical form, presented 
and extensively illustrated in [1], which has the advantage that it affords a way 
of gaining from the data being analyzed a degree of confirmation for the under- 
lying assumptions used. 


4. AN INFERENCE PROCEDURE BASED ON A Mopvutus-Ratio Statistic 


Let y; , *** Ym denote, as above, the m = N — 1] = 2? — 1 independent 
normally-distributed “contrasts” (estimates of “effects” other than the grand 
mean), having common unknown variance o’, obtained in a single replicate 
of a 2” experiment. 

Let the absolute values (moduli) |y,|, arranged in non-decreasing order, be 
denoted respectively by y/ , yi , --- y . Let a denote the integer closest to 
(.68) (m + 1), and let t,, = y{/y, . For that case of the inference problem stated 
in Section 3 which is defined by taking r = 1 (i.e., at most one of the statistics 
Yi, °** Ym has @ non-zero mean), we consider statistical inference procedures 
based on the statistic ¢,, and having the following form: 

Let k be a suitably chosen positive constant. If the observed value t,, < k, 
infer that the means of all the y,’s, 7 = 1, --- m, are zero. If t,, = k, infer that 
the true mean underlying the observed y; whose modulus is |y;| = yJ is not 
equal to zero. 

For the stated problem, such an inference procedure is recommended by the 
following considerations, in addition to plausibility and simplicity: 

1. It can be carried out in a simple, convenient, and natural way within the 
above-mentioned graphical procedure of [1]. This graphical procedure is useful 
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for drawing inferences of somewhat broader scope than jn represented in 
the inference problem stated above, including a possible graphical confirmation 
of the formal underlying assumption stated above. 

2. It is an approximate form of the optimal inference procedure for the 
stated problem, which is given in Section 7 below. 


The following Section 5.is devoted to the quantitative statistical properties 
of the inference procedure based on ¢,, . 


5. OPERATING CHARACTERISTICS OF THE MopuLus-RATIO PROCEDURE 


In order to determine the probabilities of the various possible sorts of errors 
associated with the use of the above-described procedure based on ¢,, , we con- 
“9 here the ne Reig of ¢,, soap each of the eae hypotheses: 

: E(y,) = 0 for? = 1, 


: E(y;) = Ao for just one se value j, while E(y;) = 0 for each i ¥ j; 
a A is a fixed non-zero constant. 


In order to test Hy against H, at a given significance level a, we require 
critical values k = k(m, a) such that 


1 — a = Pr {t,, < k(m, a) | Ho}. 


In Section 5.1 some approximate values of k(m, a) are given and methods of 
computation are discussed. In Section 5.2 some approximate power functions 
of such tests are given and their computation discussed. In Section 5.3 com- 
parisons are made with alternative tests based on Student’s ¢ statistic. 

5.1 CompuTaTION oF CritTIcaAL VALUES k(m, a) 


A. Case of very large m. 


When m is very large (this will be made more precise in Section C below), 
we have that yf = o with probability very near one under both Hy and H, , for 
any A; and also that ¢,, is distributed almost identically with max (ly,|, --- 
lym|)/o = yZ/o under both H, and H, . Hence 


Pr{tn <k| Ho} = Pr{ly:|/o <k for «= 1, +--+ m| Ho} 


= [&(k) — &(—k)]” = [2(8(k) — 2))", 


where 


&(k) = 


k 
= e“’”? du. 


Defining k*(m, a) by 


1 — a = [2(®(k*(m, a)) — 3)]", 


ae 1l/m 
k*(m, a) = o( + aaa"), 
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which can be computed from any full table of 6 or ~', we have that, for very 
large m, k(m, a) = k* (m, a). A few values of k* (m, a) which can be used when 
m > 60 (e.g. in 2? factorial experiments with p = 6) are given in Table 1. 


TaBLe I. Some values of k*(m, a) 


m .40 


63 (=26 — 1) 2.65 
127 (=27 — 1) 2.88 


B. Case of m = 31. 


The case m = 31 is of interest for use with 2° factorial experiments (and 
2”°* fractional replicates when p — gq = 5). To investigate the distribution of 
t, in this case, an empirical sampling study was undertaken by the Applied 
Mathematics and Statistics Laboratory of Stanford University. Table II below 
gives estimates of k(31, a) based on a sample of 2500 independent observed 
values of ts, under H, . C. Daniel has found that this empirical distribution 
of ¢,, is accurately summarized by the approximation 


Prity <k| Ho) = @(282E— 35), 


ll 
for values of kK = 2.24 (that is, for .5 S Pr {ts, S k |Ho}). 


C. Case of moderately large m. 


The effect of approximating k(m, a) by k* (m, a) is represented by the dis- 
crepancy 


D(m, a) = Prob {t, > k*(m, a) | Ho} — a. 


D(m, a) is the increase in the Type I error rate caused by using k*(m, a) instead 
of k(m, a) as a critical value. Table II gives some values of D(m, a) form = 31. 


Taste II 


k(31, «) 


REFASBaAS 
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It is seen that |D(31, a)| < .035 for a < .5. Since |D(m, a)| decreases to 0 as m 
increases for any fixed a, the values given provide a useful indication that for 
m > 31 the approximation of k(m, a) by k*(m, «) is satisfactory for most practical 
purposes; of course, for m = 31 this approximation would not be used because 
values of k(81, a) are available. 


5.2 PowER FUNCTIONS 
A. Case of very large m 
As indicated in Section 5.1 A, under H, we have for very large m 
Pr{t, <k| Ha} = Prob {ly:|/o < k, for «= 1,--- m| Hy} 
= [®(k) — &(—k)]""[@(k — A) — &(-k — A)] 
= (1 — a)" "" [(k — A) — &(—k — A)] 
= (1 — a)[@(k — A) — &(—k — A)], where k = k(m,a). 


Thus with tables of @ a simple calculation gives approximately the power 
function 


1 — 8 = 1 — B(m, a, A) = Pr{reject Hy | Ha} = Pr{t, > k(m, a) | A} 
= 1 — (1 — a)[®(k(m, a) — A) — &(—k(m, a) — A)]. 
A related operating characteristic of interest is the sensitivity function 
y = y(m,a, A) = Pr{|y;|/y2 = k(m, «) | A} 
= 1 —[&(k(m, a) — A) — ®(—k(m, a) — A)]) 
= 6(—k(m, a) — A) + &(—k(m, a) + A). 


‘y may be interpreted as the probability that the correct observation y; (for 
which E(y;) = Ao ¥ 0) will “appear as significant” (that is, that |y,|/y{, = 
k(m, «), whether or not |y;| = y,/). Clearly a/m < y(m, a, A) < 1 — B(m, a, A) 
for all A + 0 (and all m). In particular, with very large m and with A = k(m, a) 
we have 


y(m, a, A) = } + &(—2k(m,a)) > 3, and 
1 — B(m,a, A) = 1 — (1 — a)[ — &(—2k(m, a)). 
If A = k(m, a) = 1,3, we have 


y(m, a, A) se 3 and 


1 — A(m,a, A) = 3+ 5: 


Thus to obtain sensitivity y = 3 against any specified value of A = 1.3, we 
must use the critical value k = k(m, a) = A, which gives a Type I error rate of 


a = 1 — [2(@(A) — 3)”. 
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B. Case of moderately large m. 


Suppose E(y,,) = Ac ¥ Oand E(y,) = --- = E(yn-1) = 0. Let, , +++ Uns 
denote the values |y,|, --- |ym-1| arranged in non-decreasing order, let u = |yal, 
and let y{ , --+ y,, denote as usual the values u, v, , «++ Un—1 (or lys|, --* |Yml) 
arranged in nondecreasing order. Then for each sample y; , --- Ym , Since ¥.-; S 
y, Sv , we have 


U/Ve S Ym/Va < tm = Ym/YS S Yn/Ve-1 = MAX (U, Vm-1)/Ve-1 5 
and also 
Ujv, S U/Ys S tm = Ym/Yo- 

Hence we have the bounds on the sensitivity function y 
Pr{u/v, 2 k(m, a) | A} < y(m, a, A) = Pr{u/ys > k(m, a) | A} 

< Pr{u/v.-1 > k(m, a) | A}, 
and the bounds on the power function 1 — 8 
max [a, Pr{u/v, > k(m, a) | A}] < 1 — B(m, a, A) = Pr{yn/ys = k(m, a) | A} 


< Pr{u/v.-, > k(m, a) | A} + Pr{v,-./v.-1 > k(m, a)}. 
Using the facts that the density function of |y,| is 


on 
fo) = Jre"™", D0, 


and that the cumulative distribution function of |y,| is 
F(v) = aia(2) — 3], v2 0, 


we obtain (cf. for example [2], p. 369) that when m is moderately large and 


1 <i << m, v,/o has approximately a normal distribution with mean 1/c; = 
F~ (i/m) and variance 


T e«~* u(m a 2) ‘ 
2°) (m— 1)m 
From tables of @ we find, for m = 31, that 


1/e. = 1/es = r(2) = F~*(.7097) = 1.057, or cx, = .9461 


1/ea-r = 1/en r-(21) = F~’(.6774) = .989, orca, = 1.011. 


Thus C23¥22/0 has approximately mean 1 and variance .02951, and c2,v2,/0 has 
approximately mean 1 and variance .03111. 


It is useful to compare with the preceding distributions the distribution of 
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8,/0, where s, = ().{ 2?/q)' and the z,’s are independently normally distributed 
with mean zero and standard deviation o, which is for moderately large g approxi- 
mately normal with mean 1 and variance 1/2q. Thus Co2¥.3/¢ is distributed 
approximately as s,/o with q = [1/[2(.02951)]] = 17 degrees of freedom, and 
CaV2/o is distributed approximately as s,/o with q = [1/{2(.03111)]] + 16 
degrees of freedom. Similarly, for any m > 31, we have that c.v./o and ¢,_1v,_,/¢ 
are each distributed approximately as s,/o with g = [(.54)m] degrees of freedom; 
and for m > 31, .9385 < c,/ce.-, < 1, with c,/c,-, approaching 1 and also c,_, 
and c, approaching 1 as m increases. 


Hence, (yn — Ac)/c.v. has approximately Student’s ¢ distribution with q¢ 
degrees of freedom, and 


Pr{u/v, > k(m, a) | A} = Pr{u/cev, > k(m, «)/c, | A} 
= Pr{u/s, > k(m, a)/c, | A}. 


The last probability can be obtained from tables of the power function of 


Student’s ¢ test, or from the general approximate formula for that power func- 
tion given on p. 392 of [3] 


Pr{u/s, > k(m, a)/C, | A} * &((—k(m, a)/c. i A)/ Vi+ k*(m, a)/2c; q) 
+ &((—k(m, a)/c, + A)/ Vit+ k*(m, a)/2¢; q); 


which is 


= Pr{u/v, > k(m, a) | A} < . A) and 


1 — B(m, a, A). 


The quantity Pr {u/v._, > k(m, «)| A} required for upper bounds on y and 
1 — 8 may be approximated by a similar formula, or more roughly by the same 
formula found here for Pr {u/v, > k|A}. We have also Pr {v,-:/v.-1 > 
k(m, a)} = a when m is not-small. For the case m = 31, these approximations 
were used to obtain the bounds shown in Table III for the sensitivity functions 
and in Table IV for the power functions of tests at several significance levels. 


Tasis III. Approximate Bounds on 7(31, a, A) 


A 
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TaBLe IV. Approximate Bounds on 1 — (31, a, A) 


4 


4 


10 32 65 89 
18 44 ene 98 


20 21 53 .83 .97 
? 46 .80 1.00 1.00 


40 67 91 99 
0 41 17 1.00 1.00 1.00 


* These values were obtained in the empirical sampling study described above. 


It should be emphasized that these tables and the preceding discussion are 
based on the strong assumption that at most one of the y,’s has a non-zero 
mean. If the means of many y,’s are appreciably different from zero (relative 
to c), the power and sensitivity of tests based on ¢,, may be much reduced. If 
only a small number of y;’s have non-zero means, it seems likely that the power 
and sensitivity properties indicated above will tend to hold approximately, 
with appropriate reinterpretations for the various possible alternative hypotheses. 


6. COMPARISON WITH CONVENTIONAL INFERENCE PROCEDURES BASED ON 
SrupEent’s ¢ Statistic 


The standard method of analyzing a 2” (or 2”"*) factorial experiment is based 
on an underlying assumption that certain interactions are zero: H* : E(ym-a41) = 
E(ym-as2) = +++ = E(yn) = 0. In contrast, the inference procedure described 
above is based on an underlying assumption of the form: H: at most (a specified 
number) r of the m means E(y,;), i = 1, --- m, are non-zero. (The above dis- 
cussion was limited to the case r = 1.) The standard method utilizes ¢-tests 
of one or several of the hypotheses H* : E(y;) = 0,7 = 1, --- (m — d), or related 
F-tests (or related interval-estimates) to test simultaneously several of the 
hypotheses H* . In contrast, a test based on ¢,, tests (with r = 1) Hy against 
H, : E(y;) ¥ 0 for some one unknown j, j = 1, --- m. Since these problems of 
testing H* and H, are essentially different, there is no completely satisfactory 
formal basis for comparing the operating characteristics of the two types of 
procedures.. However, in applications in exploratory research experiments, it 
will sometimes be a matter of more or less subjective judgment whether H or 
H* is the more nearly appropriate underlying assumption for a given experiment; 
for cases where both H and H* seem not obviously inappropriate, the following 
comparisons of operating characteristics may be of interest. 

The quantity Pr {t,, = k(m, a)|Ho} = a represents a Type I error rate per 
experiment. If (Student’s) t-tests of g hypotheses H* : E(y;) = 0,7 = 1, --- g, 
1 <g S m — d, are made on data from a single experiment (assuming the 
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validity of the usual-corresponding underlying assumption Bly) = 
E(y.) = 0), the Type I error rate per experiment may be represented as 


Prob { at _ luel/ee = k | H%}, 


which has the value a if the critical ihe k used with the g t-statistics |y;|/s, is 
chosen suitably as a function of g, d, and a: the appropriate critical value 
k = k(g, d, a) is the upper-a point of the Studentized maximum modulus dis- 
tribution. Tables of k(g, d, .05) are available in [4]. The sensitivity function of 
this multiple ¢-test (or Studentized maximum modulus test), against an alterna- 
tive 


Hi: E(y)) = Ao #0, Ely) = +++ = Ey,) = 0, 


Pr{ly:|/s > k(g, d, a) | A} = y*(g,d,a, A), say. 
This is formally a power function of a t-test, for which (as above) we have 
*(g, d, a, A) + &((—K(g, d, a) — A)/V1 + kG, d, a)/2d) 
+ &((—K(g, d, a) + A)/V1 + kG, d, «)/2d). 


Such approximate values of y* are given in Table V for a = .05 and for several 
values of g, d, and A. Comparison of Table V with the corresponding entries of 
Table III, namely .32 < y(31, .05, 3) < .39 and .65 < (31, .05, 4) < .72, shows 


TaBLE V. Approximate Values of y*(g, d, .05, A) 


31 56 

47 77 .30 
56 86 .40 
67 .92 53 


that, at a common .05 type I error rate, the ¢,, statistic is about as sensitive in de- 
tecting one among 31 effects as: 


(a) a multiple ¢-test in detecting one among 5 possible effects, with d = 5. 
(But if 10 or more degrees of freedom can be used for error, the multiple t-test 
is more sensitive and should be used if only 5 possible effects are to be “screened.” 

(b) a multiple test in detecting one among 15 possible effects, with d = 10 
degrees of freedom. (If we assume hypothetically that 20 or more degrees of 
freedom are available for error, the corresponding multiple ¢-test would be 
preferable. However a total of g + d = 15 + 20 = 35 degrees of freedom are 
in fact not available in a 2° experiment.) 
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(c) a multiple é-test in detecting one among 20 possible effects, with d = 20 
degrees of freedom for error. 

(d) a multiple ¢-test in detecting one among 31 possible effects, with over 20 
degrees of freedom for error. 


In (c) and (d), the ¢,, test’s sensitivity is like that of multiple ¢-tests using 
more independent degrees of freedom than are actually available in a 2° experiment; 
of course this is possible just because the special form of the alternative hypo- 
thesis H, considered (r = 1) allows efficient estimation of o from the same 
degrees of freedom (y,’s) which are being screened for non-zero means. (AS 
above, it can be shown that, under H, , y/ is about as efficient an estimate of 
g as 8 with q = [.54 (m + 1)]. However, the appropriate criteria for comparing 
efficiencies here are the power or sensitivity functions of the procedures.) 

Since critical values at levels other than the .05 are not available, the above 
comparisons are restricted to this level. Higher levels would seem very useful 
with both types of procedures. 


7: A FoRMALLY OpTIMAL INFERENCE PROCEDURE 


This section is devoted to the optimal solution of the formal inference problem 
stated in Section 3 above. In principle, a most efficient inference procedure for 
this problem might be expected to have the greatest practical usefulness as well 
as theoretical interest. This seems to the writer not to be the case here, for the 
following (above mentioned) reasons: The model and assumption (e.g., r = 1) 
on which this inference problem are formally based are, in exploratory research 
situations, typically postulated somewhat tentatively as more or less highly 


plausible. Hence a solution which is optimally efficient when the exact validity 
of the underlying assumptions is postulated may not be preferred to a method 
which is somewhat less efficient in this formal sense but incorporates some 
partial check on the validity of the underlying assumptions. The latter advantage 
(which is roughly analogous to the concept of “robustness” of statistical pro- 
cedures) is afforded by the procedures of the preceding sections based on the 
t,, Statistic, as they are incorporated within the graphical methods of [1]. 

The statistical inference problem to be considered may be stated formally as 
follows: there are m normal populations, with unknown means a; and common 
unknown variance o”, from each of which a single independent observation y; 
is available, 7 = 1, --+ , m. It is known that at most r of the means may differ 
from zero (where r is given, 1 < r < m). It is required to infer on the basis of 
the observations which, if any, of the means actually differ from zero. 

(Such an inference problem could conceivably also arise as a restricted case 
of the slippage-test problem where (a) the “control state” of no slippage is one 
where the means have a common known value, (b) departures from control in 
both positive and negative directions are possible, and (c) only one observation 
is available from each population. The slippage problem is usually formulated 
as in [5] with r = 1 and without the restrictions (a) to (c). This problem could also 


be considered a restricted case of the problem of detecting “outlying observa- 
tions.’’) 
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We consider this problem first in the 
Case r = 1: The inference problem then is that of choosing, on the basis of 
one observed value of the vector y = (y; , --+ Yn), one of the m + 1 conclusions 
(hypotheses): 
dy :a; = 0, j=l,--m, 
d;:a;#0, and a;=0 foreach j #7 
foreach «= 1, +++ m. 
(For this case, the following formulation and solutions of our problem are 
analogous to those of Paulson’s treatment [5] of the slippage problem; Paulson 
considered the case of one-sided alternatives and unknown “control” mean.) 
For any inference rule § = 6(y) (i.e. any rule which, for each possible observed 
y, indicates a particular conclusion d = 4(y)), let 
Bo = Prob { d(y) = do | do true}, 
B; = Prob {&(y) = d; | d; true}, t=1,-++ ,m. 
It is natural to require 


(a) that 6(y) should be independent of the units in which the y,’s are expressed, 
and hence that 6(cy) = 4(y) for any constant c; and 
(b) that 6(y) should depend on each y; only through its absolute value |y;|. 


It follows that 5(y) should depend on y only through the statistic 


z= (a; poss » Race) — (lys/Ym Soa |Y¥m—1/Ym|) « 


Therefore we consider only inference functions of the form 6 = 4(z). Clearly, 
for any such 6, 6 is independent of oc, and 8; is determined by the value of 
la;|/o. Considerations of symmetry lead to the following definition of optimality 
of an inference rule 6(z): 


(A) Bo shall be not less than some prescribed value £3 . 
(B) For a given A > 0, the minimum of 


B; = Prob {6(z) = d; | |a,/o| = A}, for i= 1,-++-,m, 


shall be maximized. (Let 8/ denote the resulting common value of the §’s, 
which will depend on 8} and A.) 


It has been shown (the derivation is given in [6]) that the inference functions 
for this problem which are optimal in this sense are of the following form: 


iy) = d; if y? = max yj; and a/ Sv >-, 
lsism i=1 


a suitably chosen constant; and furthermore that 


iy) = dy if max i/ Su <e. 
1sism i=1 
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The suitable choice of the constant c depends on m and the chosen value 6j , c = 
c(m, 85), and is determined by the condition 


Bo = Prob { max a/ > yj < c(m, Bo) | do true} = Bo. 
lsism i=1 


Since such an inference procedure does not depend on A, it may be described as 
the uniformly best symmetric scale-invariant procedure having the specified 
Bo . Values of c(m, Bo) tabulated for a different application are available in [7]. 
(For example c(15, .95) = 0.4709 and c(15, .99) = 0.5747. The notation “k” 
in [7] corresponds to our m; we take “n” in [7] equal to 1.) With interpolation 
the tables in [7] provide values of ¢c for By = .95 and .99, and m < 120. For 
m > 120, we may use Cochran’s approximate equation (19), p. 51 in [8], which 
gives 


Bo = (1 wie gun. or c(m, Bo) = -2 log (1 oz By"). 


For a given m, a choice of 8, determines such a procedure completely, and the 
value of 1 — 8) resembles the familiar “‘significance level.” This choice can in 
principle be made with consideration of the power function 8, = 8, (A/c, Bo) 
(probability of rejecting d)) of the procedure determined by any given >. An 
upper bound on this function is given by 


Lar tx2t esx 


which can be computed by use of tables of the non-central ¢-distribution; here 
x3, *** » X2 are independent Chi-square random variables with one degree of 
freedom, and x’; , a+ is an independent non-central Chi-square random variable 
with one degree of freedom and noncentrality parameter A’. Apparently tables 
of 8; (A, Bo) are not now available. Also of interest, in principle, would be the 
sensitivity functions of such procedures, defined analogously to those given 
above for the procedures of Section 4. 

We now consider briefly the next more general case of our inference problem. 

Case r = 2: (This case will serve to illustrate the rather formidable com- 
plexity of the optimal inference rules for all cases with r = 2.) The problem 
here is that of choosing, on the basis of one observed value of the vector y = 
(Y¥1 , *** Ym), one of the 1 + m + [m(m — 1)]/2 conclusions (hypotheses): 


d,,d,,-+-d, defined as in the Case r = 1 above, and 


B,(A, Bo) < Prob — = e(m, aa, 


d;,;:a;~O and a; #0,j #1, and a =0 foreachk #7 and k ¥ j; 
i,j=1,-+-m, with ¢ <j. 


For any inference rule 6 = 6(y), let By , 6: , --- 8, be defined again as in the 
Case r = 1 above, and let 


8;; = Prob {é(y) = d,; |d;; true}, foreach 7,7 = 1,--- m, with 7 < j. 
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It is natural to restrict consideration to inference functions 6(y) which again 
satisfy the requirements (a), (b) of Case r = 1 above. It is also natural to define 
optimality of an inference function by the following conditions: 


(A) Bo shall be not less than some prescribed value £% .. 

(B’) For a given A > 0, 6; 2 8” for each 7 = 1, --- m, where 8” is a pre- 
scribed value. (Here 6’’ must be taken less than the optimal 8/ of Case r = 1 
above; or else the requirement 8; = 6” = 6{ can be met only by a procedure 
which is optimal for the Case'r = 1, and which hence never leads to an inference 
that one of the hypotheses d,; is true.) 

(C) The minimum of 


B;; = Prob {4(y) = di; | le:|/o = |e;|/o = A}, 
fort < j, 7, 7 = 1, «+: m, shall be maximized. 


It has been shown (in [6]) that the optimal inference functions 5(y) depend 
on y only through the two statistics 


u = max wil vi) 


lsism i=1 


m 4 
v = second largest wil it) ’ 
i=1 


1sisgm 


and that they have a form illustrated in Figure 1, with the d,; there having 
the same subscript as the largest y? , and the d;; there having the same subscripts 
as the two largest y;’s. 


choose d, choose qd; 


FiaureE 1 
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This procedure appears to have a form too complicated to be used in any 
anticipated application. For the operating characteristics and power or sensitivity 
functions of such a procedure, extensive and special new tables would be required. 
Furthermore, the choice of a particular case of the procedure requires a decision 
as to the relative importance (or cost) of an error when some d, is true as against 
when some d,; is true; the basis for making such a decision in quantitative 
terms is far from clear in most application situations. Similar remarks apply 
with greater force to all cases with r > 2. These complications illustrate some 
of the considerations, referred to above, which seem to indicate the preferability 
for typical exploratory research applications of procedures related to the ?,, 
statistic described in Section 4 above and in [1]. 
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Quality Control Methods for 
Several Related Variables 


J. Epwarp JACKSON 


Virginia Polytechnic Institute* 


The methods described in an earlier article devoted to control methods for two 
related variables are extended to the case of more than two related variables. The 
concept of matrix notation is introduced because of the resultant simplification in 
multivariate analysis and the original two-variable problem is restated in matrix 
form. The method of principal components is introduced both as a method of charac- 
terizing a multivariate process and as a control tool associated with control pro- 
cedures. These methods are illustrated with a numerical example from the field of 
ballistic missiles. Approximate multivariate techniques, designed to simplify the 
administration of these control programs, are also discussed. 


INTRODUCTION 


Many quality control operations in industry consist of making more than 
one type of measurement on a particular inspected item since there is more 
than one criterion for the specification of its quality. If a particular item has 
two different measurements made on it, then there exist two numbers, say x 
and y, with which the quality control personnel must determine the answer 
to the question, “Is the process in control?” The procedure most generally 
employed in the past has been to set up one control chart, say a control chart 
for averages, for each variable. Usually the process is said to be “out-of-control” 
as soon as one or the other control chart shows an out-of-control condition. For 
the discussion that follows, x and y are assumed to be normally distributed. 

In an earlier article entitled “Quality Control Methods for Two Related 
Variables” [1] it was shown that, when two variables are being controlled simul- 
taneously, the control region should take into account the relationships between 
them. The use of separate control charts for related variables, in addition to 
yielding conflicting answers, will also result in incorrect probability statements. 
For instance, the Type I error, a, is the probability that a point on the control 
chart would be out of control when the process is exactly on standard. If 95% 
control limits are being employed on a single control chart and the process is 
exactly on standard, on the average .05 of the plots on the control chart would 
be outside the limits and .95 would be inside. If there were now two control 
charts for two related variables x and y, the probability that each one would 


* Research sponsored by the Office of Naval Research, Department of the Navy: Contract 
Number: NONR-2353(01). Task order NR 042-019 with the Virginia Polytechnic Institute. 
Reproduction in whole or in part is permitted for any purpose of the United States Government. 
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plot in control when the process was on standard would still be .95 but the 
probability that both would plot in control at the same time is (.95) (.95) = .9025 
so that the true Type I error is about .10. The true control region in two-space 
is not a square or rectangle but an ellipse with all points on its perimeter having 
equal probability of occurrence. The situation is more complicated when the 
variables are correlated. In this case, the control region is still an ellipse but 
it is rotated to take into account the correlation between the variables. 

The major axis of this ellipse is called the orthogonal regression line in that it 
minimizes the sums of squares perpendicular to the line. The length of the 
semi-major axis is equal to ~/\,7'; and the length of the semi-minor axis is 
equal to ~/),7'2 where X, and ), are the roots of the equation: 

A= [(s + 8)  V(s; + 8)” — Asis, — (8,2)°]]/2 
and 
_ AN — IF. 
: ~ N-2 
(where F hasn, =,2 and n. = N — 2 degrees of (freedom), s2 , s} and s,, having 
been obtained from a base period sample of size N. 

An example was given in this earlier article where such an ellipse was con- 
structed based on the variability of seventy-five observations from a chemical 
process in which analytical determinations were made each time on two con- 
stituents (x = hydroquinone and y = Elon; both in terms of [grams/liter] - 100). 
From these sample measurements, the following quantities were obtained: 


& = 264.32 g = 471.48 
s; = 102.65 s; = 107.96 3,, = 68.87. 


From these quantities, a control ellipse was constructed using the sample average 
as the center. It was shown that this method of control was identical to that of 
Hotelling’s T°-charts [2] in that all points inside the ellipse would be “‘in-control” 
on the T?-chart, all points outside of the ellipse would be outside the 7?-limit 
and all points on the perimeter of the ellipse would also be located on the T’-limit 
line. The equation for the ellipse is: 


qe — tee | — 4" , y- i" _ 28a — Bly — ). 


T? 


ae 2 2 
$8, — Szy 8, Sy 8,8, 


In some control situations, arbitrarily imposed standards are used in place of 
the sample means. However, in this article the standard values will be assumed 
to be the sample means. 

This procedure amounts essentially to transforming the original two correl- 
ated variables into two new variables that are independent. The method is 
known as principal axis rotation or the method of principal components [3]. The 
axes of the ellipse may be thought of as vectors in two-space and in fact are 
referred to as characteristic vectors in that they “characterize” the rotation 
required of this two-dimensional problem to obtain independent variates. 
(Engineers may find the synonym eigenvectors more familiar.) Similarly, the 
quantities \, and \, are called characteristic roots. (Latent roots and eigenvalues 
are synonymous terms.) For this chemical example, the vector components 
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corresponding to the major axis were v,, = 9.151899 and v,, = 9.511514; the 
components for the minor axis were v2, = — 4.346576 and v2, = 4.182239. The 
characteristic roots were \, = 174.22616 and A, = 36.38384. The remainder 
of the article consisted of a full description of the numerical example including 
some computing techniques and checks. It is the purpose of this present article 
to illustrate how these characteristic vectors may themselves be employed both 


to describe the nature of-a multivariate situation and to control the process 
that it represents. 


Tue Ust or Matrix ALGEBRA 


In the appropriate situations where matrix notation may be used, it has 
proven a boon to authors, printers and readers alike in the resultant simplifi- 
cation of mathematical presentations. It is perhaps appropriate at this point 
to restate some of the material of the previous article in matrix terminology 


in order to provide a more uniform approach for the remainder of the present 
article. 


Matrices will be designated by boldface upper case characters. Row vectors 
will be designated by lower case boldface characters. A prime denotes matrix 
transposition (i.e., A’) and A™* refers to the inverse of A. TrA refers to the 
trace of the matrix A which is the sum of its diagonal elements. 

Let the sample covariance matrix be 


S = |" ‘| 
en & 
The characteristic roots \, and 2 can then be obtained by solving the de- 
terminantal equation: i 


2 ® 
is —al| = F oe it 
8, 8 —X 


Hb} ky 
0 r». LO 1 


the identity matrix, and let v, = [v,.0,,] be the first characteristic vector and 
Vz, = [v2,Vey] be the second characteristic vector so that 


“arr 
v2 Vo, Vo 
In a similar manner, the computational checks may be rewritten: 
l.twA=#S, ie. AtAW=K+8, 
2. |A] = |S], ive. AiA2 = 838; — (s,,)’, 
3. VV’ =A, i.e. vy = Ay , V2 = Ao, 
v3 = 0, 
4.VV=S, ie. wt, +03,= 82, ete. 
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Finally, letting x = [(x — #) (y — 9)], 


qt = ze [58 Wo _ male — A — 9) - sy. 


88, — (82y) 
At this point, it would be best to depart a bit from the previous notation 


by letting 7, = x — Zand z, = y — 9. (ie., x = [2,22] = [x — Z, y — 9]). This 
will facilitate the extension of these methods to more general cases. 


Tue Principat Axis Rotation 


As stated in the introduction, in order to determine the major and minor 
axes of the control ellipse, one performs a principal axis rotation on the original 
variables x, and 2x, transforming them into new variables z, and z, that are 
uncorrelated. This transformation is merely: 


[2:20] = [x,2,.])V’ or z= xV’. 


The covariance matrix of the z’s is obtained by pre- and post-multiplying 
the covariance matrix of the z’s by the matrix of the transformation V (viz. 


VSV’). This quantity is equal to 
Dae 
9 


For the example mentioned in the earlier article, where x, = hydroquinone 
and z, = Elon, VSV’ = 


9.151899 or ie pei 


—4.346576 4.182239JL 68.87 107.96L9.511514 4.182239. 


ais 0 
0 1828.78 


The new variable z, will have zero mean and variance }? . z, will also have zero 
mean and variance \3 . The covariance between z, and z, is zero, hence z, and 2, 
are uncorrelated. The original variables x, and x, are assumed to be bivariate- 
normally distributed and if this is true, the new variables z, and 2, will also 
be normally distributed since this is a linear transformation [4]. 

To use this as a control tool, it is sometimes convenient to divide each element 
Of the characteristic vectors by its corresponding characteristic root. That is, 
let w, = v,/A, and w, = v2/Az so that the transformation on the original variables 
is now y = zA~* = xW’. Then the covariance matrix of the new variables, the 
y’s, is: WSW’ = I, the identity matrix. (see Appendix A). In other words, the 
new variables, y, and y2 , in addition to being independent and being bivariate- 
normally distributed with zero means, now also have unit variances. For example, 


W,= [w1:%,2] = [.0526 0546], 
WwW, = [we: Wee] ed [—.1195 1149], 
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and 


wsw’ “| 0526 eae ce ee ao a ' 00 | 
—.1195 .1149]L 68.87 107.96]L.0546 = .1149 00 1.00 


The sum of squares of p independent, normally distributed, variables with zero 
means and unit variances has the chi-square distribution with p degrees of 
freedom (in this case, 2 degrees of freedom) [5]. For small samples such as this 
(N = 75), when the variances are estimated from the sample, x’ must be replaced 
with 7”. The sum of squares of individual pairs of y’s is equal to T? (i.e., T? = 
yy’, see Appendix B). 

From the same previous example: let a sample observation be 2, (hydro- 
quinone) = 280.32 and z, (Elon) = 463.08. Then, since Z, = 264.32 and Z, = 
471.48, x = [x, — #, 2, — #,] = [280.32 — 264.32 463.08 — 471.48] = [16.00 — 
8.40]. 


T? = xS"x’ = [16.00 — $40] AN oo pa = 8.43. 
— .0109 .0162jL—8.40 
If transformed variables are used: 


0526 Brn) = [.3830 — 2.8772]. 


.0546 1149 
That is, y, = .3830 and y, = — 2.8772. 


y = xW’ = [16.00 — $40] 


yokes = (.3830)? + (—2.8772)° 


— 2.877 = 8.43 


which demonstrates the equality of the two methods. Since the 95 per cent 7” 
limit for two variables based on a sample size of 75 is 6.325, the process may be 
considered to be out of control. The point would also plot outside the ellipse 
shown in the previous article. 


7? = yy’ = [.3830 — 2772 


ALTERNATIVE METHODS OF CONTROL 
It can now be seen that three alternative methods of control are available: 


1. Plotting the original observation (x; , 72) on the control ellipse. 
2. Using a T” chart for the quantity xS~'x’. 
3. Using a T° chart for the quantity yy’. 


Methods (2) and (3) will, of course, yield identical results on the 7” chart. 
Even though Method (3) must be carried out in two parts (i.e., y = xW’ and 
T, = yy’), the computations are about the same as for Method (2). All three 
methods will give the same answer to the question “Is the process in control?” 
Which of these three is the best to use? 

One of the conclusions of the previous article was that Method (1) would 
indicate, graphically, the nature of out of control conditions while Method (2) 
land now Method (8) also] would preserve the time scale and numerically 
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summarize the condition of the process by one number. Method (3) can in 
addition, however, be used to describe the nature of the process as follows: 

The coefficients of W associated with y, in this example are both positive 
and of roughly the same magnitude (.0526 and .0546). Therefore, y,; = 2,w,, + 
LW. = .0526z, + .05462, might be zonsidered as a measure of average process 
variability (i.e., y, will fluctuate in the same way as the sum, and hence, the 
average, of x, and x,). On the other hand, the coefficients of W associated with 
Y2 , although they too are of the same order of magnitude, have opposite signs 
(—.1195 and .1149). Therefore, y, = —.11952, + .11492, might be considered 
a measure of the disagreement between the two variables x, and x, regarding 
the nature of the process. It may be remembered that the sum of the character- 
istic roots A, and A, equals the sum of the variances s? and s> . Thus, s? + s° is 
the trace of this matrix S$ (written érS) and can be thought of as a rough measure 
of the total variability of the system. The ratio of each of the characteristic 
roots to trS will give an estimate of the proportion of the total variation explained 
by each characteristic vector. In this case: — 


d,/trS = 174.23/210.61 
d./trS = 36.38/210.61 = .173. 


This implies that about 83 per cent of the totalivariability is due to the “process” 
or things common to both analyses and the remaining 17 per cent is due to the 
disagreement between analyses. 

As long as the T’-chart is in control, usual operation may be continued. As 
soon as the T’-chart is out of control, reference should be made to the correspond- 
ing values of y, and y, to determine the nature of the out-of-control condition. 
In the example above (remembering that both y, and y2 have zero means, unit 
variances, and hence unit standard deviations), y, = .3830 so the process itself 
would seem to be in control but y, = — 2.8772, which would be outside two- 
sigma limits for this variable, and would indicate that there is a disagreement 
between the hydroquinone and Elon analyses. Referring back to the original 
data, it can be seen that the hydroquinone is above standard and the Elon is 
below standard. Since these two variables have high positive correlation, a 
difference as large as this cannot be expected by chance alone and hence this is 
the cause for the T’-chart going out of control. Method (3), then, is a combi- 
nation of Methods (1) and (2) in that not only can it give one answer to the 
question, “Is the process in control?” but it can also indicate the nature of out- 
of-control conditions when they do exist. It must be remembered that the 
coefficients of W are only estimates of the true coefficients for these rotated 
axes and, particularly for small samples, these estimates may be misleading. 

In summary then, for the two-dimensional case, T’-control can be maintained 
either by use of a control ellipse or by use of a 7-chart with or without supple- 
mentation by control charts for the transformed variates. Method (2) would 
seem to have no advantage over the other two methods and the choice between 
the remaining two is strictly a matter of preference with regard to the individua! 
control situation. The control ellipse is certainly. much easier - to: maintain, 





QUALITY CONTROL METHODS FOR SEVERAL RELATED VARIABLES 365 


particularly in the process area, while Method (3) maintains the time scale 
and may be useful for long-range control work. 

For best results, the practice of performing the operations involved in Method 
(3) on the sample covariance matrix should be done only when all of the variables 
are in the same units and have roughly the same variances. If this is not the 
case, the operations should probably be performed on the correlation matrix. 
The variables x, and 2, will then have to be divided by their respective standard 


deviations before going through the 7” calculations and the transformation 
from the z’s to the y’s. 


EXTENSION TO THREE OR MorRE VARIABLES 


For more than two variables, there is no longer the choice in methods since 
plotting points in three-space is rather difficult and for four-space and higher, 
impossible. It is in these cases that the use of transformed variates really 
becomes helpful. Furthermore, the need for multivariate control becomes 
increasingly necessary with increasing numbers of variables. It may be recalled 
that for the two-dimensional case, where +20 (95.4 per cent) limits were em- 
ployed for each variable, the probability of a sample point being in control for 
both variables when the true means are both on standard is (.954) (.954) = .910. 
For three variables, this would be (.954)* = .868; for four variables, (.954)* = 
.828; and so on. For, say, a ten-variable problem, the probability of a point 
being in control for all ten variables when all ten true means were on standard 
is (.954)"° = .624. In other words, in a ten-variable problem (which is not 
uncommon in industry), the use of ten individual two-sigma control charts 
will result in an apparent out-of-control condition calling for corrective action, 
roughly, one time in three when actually the process is on standard in all respects. 
This sort of situation could result in such strained relations between production, 
quality control, and perhaps maintenance that the whole control program 
might be in jeopardy. Increasing the width of the limits to, say, +30 will reduce 
this error but at the price of increasing the Type II error. As in the two-di- 
mensional case, the use of a multivariate control procedure will allow the Type 
I error to be fixed for the system as a whole. However, the warning concerning 
the sampling variation of the w;,;’s, mentioned in the two-variable case, should 
be heeded even more when several variables are being considered. There is also 
no guarantee that any physical significance can be attached to some or all of 
the sample characteristic vectors. Hence good judgment is required in the 
interpretation of these vectors and their application to process control. However, 
subject to these criticisms, the method of principal components can be very 


useful in multivariate quality control work and has already proven to be so in a 
number of cases. 


Usr oF MULTIVARIATE QuALITY CONTROL 
IN BAuuistic MissiLeE TESTING 


For an illustrative example, a four-variable problem dealing with ballistic 
missiles will be discussed. One of the important properties of the performance 
of a ballistic missile is the impulse that it produces during firing. An estimate of 
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this quantity is usually obtained by a method known as static testing in which 
a loaded rocket is fastened down to prevent its flight during firing. Thrust 
gages are fastened to the head of the rocket to record the thrust exhibited. 
The total amount of thrust produced during the time in which the powder 
inside the rocket is burning is the total impulse of that particular missile. Since 
missiles such as the Nike and Honest John are relatively expensive, it is desirable 
to obtain as reliable information from each test firing as possible. For this 


x 
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Figure 1—Ballistic missile test. Control chart for original variables. 95% control limits. 
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reason, two identical thrust gages are attached to the head of the rocket. Each 
gage is connected with a separate system of recording apparatus. Each of these 
recording systems consists of (1) an electronic integrator which automatically 
determines the total impulse and (2) an oscilloscope and camera which produce 
a photographic record of the rocket’s thrust as a function of time. This photo- 
graphic record can then be planimetered to obtain the total impulse. There 
are, then, from a single firing, four impulse measurements: an integrator and 
planimeter measurement corresponding to each of the two thrust gages. Assuming 
no systematic biases on the part of any of the components, the best estimate 
of the total impulse would be the average of all four measurements. However, 
by treatinz this as a four-variable problem, considerably more information can 
be obtained regarding the reliability of the measurement system. Analogous 
to the two-variable example where two related variables were transformed 
into two independent variables, the present four-variable system can, by a 
similar principal axis transformation, be transformed into a system with four 
independent variables. 
Let the original variables be: 


= Gage No. 1—integrator reading, 
= Gage No. 1—planimeter measurement, 
= Gage No. 2—integrator reading, 
= Gage No. 2—planimeter measurement. 


Some typical control data for 15 rounds fired since the base period for this 
study are shown in Figure 1. 


Over a base period of approximately 40 test firings, the covariance matrix 
(coded somewhat for security reasons) was as follows: 


102.74 88.67 67.04 54.06 
88.67 142.74 86.56 80.03 
67.04 86.56 84.57 69.42 
54.06 80.03 69.42 99.06 


Ss = 


For more than two variables, the characteristic vectors can be obtained more 
easily by numerical methods than by direct solution. There exist a great number 
of procedures for use in problems such as these. The choice of procedure will 
depend on the size of the matrix and the type of computing equipment available. 
In the present example, the method of iteration [6] was employed. As a result of 
the principal axis transformation, the new independent variables are: 


Yi = 02562, + .03322, + .02517, + .02452, , 
Y2 = —.0897x, — .02582, + .02007;, + .108iz, , 
yz; = —.10552, + .1403z7, — .03107; — .0483z2, , 
Yo — .06427, — .03612, + .21277, — .1015z, . 
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The characteristic roots associated with these equations are: 


A, = 335.35, A. = 48.04, As = 29.33, Ay = 16.42. 


As in the two-dimensional case, the sum of the roots add up to the trace of the 
covariance matrix and, again, the ratio of each root to the total can be taken 
as an estimate of the proportion of the total variability explained by the corre- 
sponding characteristic vector. In this example, the first vector explains 78.2 
per cent; the second, 11.2 per cent; the third, 6.8 per cent; and the last, 3.8 
per cent. 


The physical ‘characteristics’ associated with each of these vectors may be 
thought of as follows: 


(1) The coefficients of y, are all positive and of roughly the same magnitude. 
y, will bear a close relationship with the average of 2, , x. , x; , and x, and can 
be thought of as a measure of product variability. In many industrial appli- 
cations this might represent not only product variability but testing variability 
that was common to all the original variables. In this example, for instance, if 
it were possible for atmospheric conditions to influence the gages, this might 
then affect all four measurements in the same way. 

(2) The coefficients of y, are positive for Gage No. 1 and negative for Gage 
No. 2. Therefore, this transformed variable may be thought of as a measure 
of the difference between gages. 

(3) The coefficients of y; are all negative except x. . However, the coefficients 
of x, and 2x, are larger than the others and hence y; may be thought of as a 
measure of the difference between integrator and planimeter results on Gage 
No. 1. 

(4) Similarly, y, may be thought of as a measure of the difference between 
integrator and planimeter results on Gage No. 2. 


It may appear to the reader that some of these definitions are quite arbitrary, 
particularly with regard to the exclusion of some of the original variables, but 
by use of a technique described in Appendix C, it will be shown that these 
diagnoses do have some justification. These results would indicate that y. , y3 , 
and y, represent various kinds of measurement variability and that at least 21.8 
per cent of the total variability is due to measurement variability, while the 
remaining 78.2 per cent is product variability plus any testing variability that 
affects all four measurements the same way. 

As before, the sum of squares of the four y’s (yy’) for each round fired can be 
plotted directly on a T?-chart as shown for the same 15 rounds in Figure 2, and 
as long as this chart remains in control, it may be concluded that the variation 
in thrust measurements is not significantly different from that of the base period. 
The four y’s may also be plotted on supplemental control charts (Figure 3) 
and can be used both to diagnose trouble when the 7°-chart indicates an out-ot- 
control situation and possibly to indicate trends which would eventually lead 
to the T?-chart going out of control. 

It should be emphasized that these four control charts should be used only 
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i 
Figure 2—Ballistie missile test. T? control chart. 95% control limit. 


in this connection, however. They should not be used in substitution for the 
T’-chart since they suffer from some of the same problems in connection with 
the Type I error as the control charts on the original variables. 

From time to time after this control procedure has been established, the new 
control data should be investigated to make sure that the relationships between 
the original variables have not changed perceptably. If the relationships have 
changed, the control procedure should be modified accordingly. 


APPROXIMATION PROCEDURES 


Although the computations involved in obtaining the y’s and T” are fairly 
simple, they may, nevertheless, be a little too much for a control man on the 
scene to handle, particularly if he must analyze a large volume of control data. 
One solution for this is to simplify the calculations so that they can be done 
without the aid of a calculator. Simplification usually means loss of precision 
but, if care is taken, the results may still be good enough that they are adequate 
for use in control work. The methods employed in this section will be described 
more fully in Appendix C. 

Earlier, in a discussion of the characteristics of the four vectors associated 
with the ballistic missile problem, it was stated that y, , the transformed variable 
associated with the first characteristic vector, would bear a close relationship 
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Fiaurs 3—Ballistic missile test. Control chart for transformed variables. 95% control limits. 


to the sum or average of the original variables. One might be tempted, for on- 
the-spot control work, to use either the sum or average of x, , 22 , %; and 2% in 
place of y, . Similarly for y. , one might want to use the difference between 
(a, + 22) and (x3 + 24) etc. In this way, a new set of transformations may be 
derived which are similar to the characteristic vectors but are much easier to 
use in control work. In this case, the transformations are: 
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WM=U1+etutry 
Us, = (a, + 22) — (ts + 24) 
Us; = X, — Le 
U, = %3 — M% 


Just picking a simple transformation is not enough. It must pass three re- 
quirements: 


(1) Orthogonality. For this example, the matrix of coefficients of these new 
transformations (call it W*) is: 


1 1 
1 1 
1 =] 
0 0 i =~] 


so that the transformation may be written u = xW*’. One requirement of the 
simplified transformation is that the sum of cross-products of the coefficients 
of any row of W* with any other row should be zero. [i.e., Row 1 X Row 3: 
(1)(1) + (1)(—1) + (1)(0) + (1)(0) = 0.] (This requirement is only one of the 
conditions which ah orthogonal matrix must fulfill, but it is sufficient for control 
purposes.) This requirement is automatically fulfilled by the V and W matrices 
by the manner in which they are obtained. 

(2) Correlation with Characteristic Vectors. Not only should the approximation 
u; bear a strong resemblance to the corresponding y, , but it should be highly 
correlated with it as well. In this example, over the base period, the correlation 
between y, and u, was .998, essentially perfect. The correlation between y, and 
u, was —.90; between y; and us , —.87; and between y, and u, , .79. (Negative 
correlations are just as good as positive correlations in these instances since 
the choice of signs in W* is arbitrary.) Generally speaking, it becomes increasingly 
difficult to get simple approximations which have very high correlation with the 
characteristic vectors as the sizes of the characteristic roots get smaller. The 
fact that the coefficients of W are themselves sample estimates will of course 
affect these correlations and, for small N these correlations may not be all that 
one might desire. 

(3) Independence. Since the y’s were independent of each other, the u’s should 
be independent also. In general, this is not possible and the goal, usually, is to 
keep the correlations as small as possible. In this example (see Appendix C) 
only one correlation was significant, that between u, and u, , the correlation 
coefficient being .39. Again, for small N these correlations may not behave as 
one might wish. 


Undoubtedly, some of the correlations of the approximations with the 
characteristic vectors could be increased and some of the correlations of the 
approximations with each other could be decreased by choosing some other 
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transformations but the present.set was felt to have the best combination of 
desirable correlation coefficients and simplicity of calculation for use in this 
control program. In addition, as stated earlier, the coefficients of W are only 
sample estimates and if the study is based on a small sample, the approximations 
may conceivably be fully as valid as the sample characteristic vectors themselves. 

The use of linear combinations of variables is frequently employed in industry 
but usually without the preliminary analysis of the data described in this article 
and, therefore, these three requirements are seldom fulfilled, and in fact, no 
information about the conformity of a system to requirements (2) and (3) 
would be available. If the linear combinations themselves were highly correlated, 
then the control program is little or no better than with the original variables. 

The standard deviations for the w’s, as determined in the appendix, are as 
follows: 


8, = 36.3, &,, = 13.0, &,, = 8.3, 8, = 6.7. 


There is still the problem of approximating 7°. It may be recalled that the 
sums of squares of the four y’s (i.e., yy’) has a 7°-distribution. If the w’s were 
independent of each other, then 

= E - a" 


t=1 Sus 


would also have a 7°-distribution. In this example, this quantity would equal 


(uw ee il;)” ee (us oe tly)” s. (u; — tls)” ie (Us — tis)” 
1320.67 169.91 68.14 44.79 

If the individual u,’s do not assume many different values, the quantities of 
(u; — @,)’/s?, may be tabulated and the calculation of the approximation of 
T’, call it T?, reduces to simple addition. The more the u’s are correlated, the 
poorer this approximation will be. If the w’s are positively correlated, T% will, 
in general, be less than 7”. If the u’s are negatively correlated, T? will, in general, 
be greater than 7’. Figure 4 shows the relationship between T’ and T? for the 
same data as shown in Figures 1-3. In this particular case, the relationship 
is a good one, indicating that the u’s are good approximations for the y’s. There 
is a slight tendency for T? to understate 7’, and, in fact, the one round that 
was just outside the T” control limit is just inside the T? limit. This is due to 
the fact that most of the correlation coefficients between the u’s are positive. 

It should be emphasized that these approximations and in particular, the 
approximation to 7”, are designed for use at the shop level where it is felt that 
an approximate multivariate control tool is better than none at all. Where 
the proper facilities are available (and this amounts only to a desk calculator 
and someone familiar with its operation for these calculations), the regular 
T’-calculations should be used. If the approximation techniques are employed, 
the regular T’-calculations should still be carried on for long-range control 
analysis. 

It has been suggested above that the relationships among the original variables 
should be checked from time to time. Similarly, the relationships between the 
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Ficure 4—Ballistic missile test. Relationship between T? and T3. 95% control limits. 


transformed variables and their corresponding approximations should also be 
checked. 


EXTENSION TO MorE VARIABLES 


It has now been demonstrated that the method of Principal Components can 
be used both to analyze a multivariate process (similar, say, to components of 
variance in analysis of variance) and to control the process (an extension of the 
common control chart). 

Similar methods may be worked out for any number of variables. As the 
number of variables increases, so, of course, does the amount of computation. 
However, the type of concern usually faced with a control problem consisting 
of a large number of related variables usually has access to some sort of com- 
puting equipment, and even the most modest computational installation will 
make quick work out of most of these calculations. (The computational re- 
quirements for the ballistic missile example were small enough that it took less 
than a day on an ordinary desk calculator.) 

As the number of variables increases, it is also possible that redundancy may 
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begin to appear. That is, there are not as many actual variables being accounted 
for as there are measurements being made. For example, in a ten-variable 
problem, it may turn out that the first five characteristic vectors explain over 
99 per cent of the total variability and this number is all that is necessary for 
control purposes. Although the theory associated with this type of situation 
has not been completely developed, those interested°-may wish to refer to [7] 
for a discussion of the motivation for such a situation and an approximate 
method for its solution. 

As in the univariate case, considerable more insight into the nature of the 
process can be obtained by using rational subgroups wherever possible. For a 
discussion of this, see [2], [7]. 


ApPpENDIXx A 


Proof that WSW’' = I 


Given a set of x-variates with covariance matrix §, ; if linear transformation 
of the form y = xW’ is applied, the covariance matrix of the resultant y’s is 
S, = WS,W’. By definition, W = A7'V. Then: 


WSW’ = A'VSV’A™ 
= A'VV’'VV’A” Since S = V/V 
= A“AAA™ Since A = VV’ 
=I 
Therefore, the new variates of the transformation y = xW’ are independent 


and have unit variances. 


AppEenpIx B 
Proof that T? = yy’ 
By definition, 
y = xw’ 
From Appendix A, 
WswW’ =I, S=W'Iw'=W'w", s"*=wWww 


Therefore 


T* = xS"'s’ = xW’Wr’ = yy’. 


APPENDIX C 
Approximation Procedures 


To facilitate printing operations, the number of digits in this numerical 
example has been considerably reduced. Hence, the example may not always 
check out exactly. 
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As stated in Appendix A, if a set of variables is transformed by a set of linear 
transformations, the covariance matrix of the new variables is the result of the 
original covariance matrix being pre-and post-multiplied by the matrix of the 
transformation. In the ballistic missile problem, the principal-axis transformation 
y = xW’ is represented by the matrix: 


0256 .0332 0251 0245 
— .0897 —.0258 .0200 1081 
— .1055 1403 —.0310 —.0483 ) 
— .0642 —.0361 2127 —.1015 


W= 


Since the covariance matrix of the original variables was: 


102.74 88.67 67.04 54.06 


88.67 142.74 86.56 80.03 
’ 
67.04 86.56 84.57 69.42 
L 54.06 80.03 69.42 99.06 


the covariance matrix of the transformed variables is: 


1.00 00 .00  .00 
00 1.00 00 .00 
00 .00 1.00 .00 
00 00 .00 1.00) 


WSW’ = 


the identity matrix I. This shows that the y’s are independent and all have unit 
variances as they should by definition. 

To determine, for example, the ‘relationship between y, and its approximation 
U, , form the matrix: 


> Wer Was Bt - | W2 i “~ .—.0258 .0200 sae 
1 1 -1 -1 


wh Wh Wh ws wi 


Forming the product 


[=2-|s[~ + 1 a 
| wi wi —11.70 169.91 
it can be seen that the variance of y, is 1, the variance of u, is 169.91 and their 


covariance is — 11.70; hence, the correlation between y2 and u, is —.90. Extend- 
ing this to all four vectors and their approximations, form 


(8) 
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1 1 
1 —1 
a ey ® 
ae 
This quantity is equal to 
wsw*’ 
Ww* SWw* , 
0 36.3 


Oo 5.2 44.8): 


WSW’ is the covariance matrix of the y’s; W*SW*’ is the covariance matrix of 
the u’s; WSW*’ or W*SW’ represent the covariances between the y’s and the 
u’s. ; 


' 
! 
| 
| 
1 
! 
! 
' 
| 
! 
| 
t 
| 
| 
| 
| 
! 
! 
! 
| 
| 
! 
' 


The corresponding correlation matrix is: 
1.00 00 .00- :00 | 1,00 

1.00 .00 00: 

00 1 


12 1.00 


The upper-left quadrant of this matrix should equal the identity matrix I 
by definition. The upper-right (or lower-left) quadrant represents the correlations 
between the y’s and w’s. For these approximations to be satisfactory, the diagonal 
elements of this quadrant should be close to + 1 and the off-diagonal elements 
should be close to zero. The lower-right quadrant represents the correlations 
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among the approximations themselves and these should be as small as possible. 
In this manner, any linear combination of the original variables can be analyzed 
with regard to its usefulness in a control system. In this present example, the 
diagonal elements of the upper-right quadrant are satisfactorily close to either 
plus or minus one but some of the off-diagonals are a little large, particularly 
between y, and u, . In the lower-right quadrant, none of the off-diagonal elements 
are significant except between u. and wu, . As has already been stated, better 
approximations could have been found but for a sample size of forty this set 
was felt to be a good compromise as far as matching the set of characteristic 
vectors and being easy to administrate. The over-all results are satisfactory as 
can be seen in Figure 4. 

The approximations in this example were obtained primarily by trial and 
error. Better techniques are being introduced into industrial situations and 
with the increased use of factor-analysis techniques, the methods of obtaining 
these approximations should become more straightfoward. 
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of Row-Column Interaction 
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A serious limitation in the use of latin squares is the confounding of the main 
effect of each. factor with the interaction of the remaining two factors. In some 
cases, the interaction of rows and columns can be expressed as a multiplicative 
term of assigned factors associated with rows and columns. The analysis of such 
designs is presented in detail. An example is discussed. The method is extended to 
the case where the treatment effects are linearly related to a given concomitant 
variable. A brief discussion is given of the relation of the proposed method with 
other tests for non-additivity. 


INTRODUCTION 


The usefulness of latin square designs is severely limited by the confounding 
of interactions with main effects. As an example, consider a study of the abrasion 
resistance of leather in terms of its dependence on the relative humidity main- 
tained during the conditioning of the hide prior to test [1]. A six by six latin 
square was used in this study, as shown in Table I. The rows and columns of 
the latin square. are the rows and columns of the rectangular pattern according 
to which the specimens were cut from the hide. The treatments, represented 
by the letters A through F, are the six values of relative humidity. The com- 
parison of the 6 averages corresponding respectively to the sets of specimens 
marked A, B, C, D, E, and F is entirely valid, even when non-linear trends 
exist in the abrasion resistance of the leather, both in the horizontal and vertical 
direction, provided that these trends do not interact. The meaning of “no 
interaction” is that the abrasion resistance of the untreated leather at each 
point of the hide can be represented as the sum of two terms, the one depend- 
ing only on the row-trend and the other only on the column trend. Expressed 
in a different but equivalent way, the absence of interaction requires that, for 
example, the difference in abrasion resistance (of the untreated leather) between 
specimens 2 and 3 is the same as between specimens 8 and 9, 14 and 15, 20 and 
21, 26 and 27, and 32 and 33. In practice, one is not always justified in making 
this rather restrictive assumption. 

The purpose of this paper is to describe a valid method of analysis of latin 
square data in the presence of row-column interaction, provided that this inter- 
action can be represented by a simple multiplicative term. 
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TABLE 1 

Design for Abrasion Study' 

3 4 5 

9 10 11 
15 18 
21 24 
27 30 


33 . 36 


: Menta in upper left of cell identifies the specimen; letter in center of cell represents 
level of relative humidity as follows: 


A=2% B=37% C=530% D=62% E=75% F = 87% 


THE MopEL 


Consider a latin square of side n, with a single observation per cell. Let y;;, 
represent the measured value of some given property for a specimen corre- 
sponding to row 7 and column j, and subjected to treatment k. 

We assume that the following model applies: 


Yu=ehtot+7; + BAB; + O% + €:i ~ 


In this equation ¢;;, is a random error. The symbols yu, p; , y; and 6, correspond 
to main effects and represent, respectively, the general mean, the ith row effect, 
the jth column effect and the kth treatment effect. The term 8A,B, is included to 
represent a row-column interaction of a simple multiplication type, by assigning 
to each row a quantitative variable A; and to each- column a quantitative 
variable B; . The values A; and B; , subject only to. the restriction below, are 
chosen either on the basis of a physical knowledge of the system under study, 
or on some plausible assumption regarding this system. The coefficient 8 is an 
unknown parameter. 

' The following assumptions are made: 


1. Dea= PA = Dn = VB = Da =0 (2) 
2. The ¢;;, are normally and independently distributed, with mean zero and 
a common variance o”, not necessarily known. 

Our problem is to obtain estimates for p; , y; , 9 , 8, and o’,:and to make 
appropriate tests of significance involving these quantities. 

It is seen that the model corresponds to a simple covariance situation. The 
formulas given below can, therefore, be obtained by applying the usual methods 
of covariance analysis [2]. 
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ANALYSIS 
1. Estimation of parameters 
The following notation will be used: 


Rk; = by Yiix = sum of all observations in ith row 
C; = x Yii, = Sum of all observations in jth column 


sum of all observations corresponding to kth treatment 
Dis = A,B; 
P, = sum of all p,; corrésponding to cells occupied by the kth treatment 
G = - R; = a C; = x T, = sum of all observations. 


Denoting a least-squares estimate by a caret (-), the following formulas 
are readily derived by covariance analysis: 


(3) 
1, +++ ,n) (4) 


j= 1,+-+ ,n) (5) 


ms Les YiikPis — x TP, 


LP LP " 


ioe n) (7) 


2. Analysis of Variance 


The estimate f, f; and 7; are orthogonal to # and the 6, and their sums of 
squares are therefore calculated by the usual formulas for two-way classification. 
On the other hand, # and the 6, are not orthogonal to each other. Therefore 
the sums of squares calculated for these quantities are not additive and depend 
on the order in which the estimates are computed. Table 2 contains the formulas 
required in the analysis of variance, and can be used for tests of significance, 
discussed in the following section. 


3. Standard Errors 
The standard error of # is given by the following formula: 
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TABLE 2 
Analysis of Variance 


Source of 
Variation Sum of Squares 


_— SS. = D0 Livin 


ce 


n? 


LR 


General Mean SS, = 


Rows SS, 


n 
XC 


Columns SS, 


n 
oT 
Treatments, k 
ignoring n 
row-column 


interaction 7 

(n » X YixPis — » T.P,) 

Row-Column = tS SS 
Interaction, n(n x dpi fa - a P 1) 
corrected for . 
treatments 


Error n? = SS, — (SS, + SS; + SS, + SS; + SS) 


(dh x YsixDis)” 
Row-Column SS, = 


Interaction, . Lis 
ignoring 
treatments 

Treatments, SS, = SS; + SS, — SS; 
corrected for 
row-column 
interaction 


The standard error of any linear combination of the treatment effects, 
L = >>. w6, , is given by the following formula: 


LG (20 Ps)" ; 
"la tae DS Le EPO] - 


In particular, the standard error of the difference between treatments k and 
k’ is: 


ee ee (P, — Pi) (10) 


no ou — LP) E75 
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Sample estimates of these standard errors are obtained by substituting for 
o, the square root of the mean square for “Error” (Table 2). 


4. Tests of Significance 


It is advisable to begin the analysis by testing the significance of the inter- 
action term. The appropriate F test, with 1 and (n® — 3n + 1) degrees of freedom, 
is: 


SS, 
- BS /@F — ia +1) sais 


Fs 


If this F is found to be significant, there is evidence that 8 is not zero. Equa- 
tion (6) then yields an unbiased estimate for this quantity. If the interaction 
term is significant, the test of significance for treatments should be carried out 
on the corrected treatment sum of squares. 


F treatments = 7 Hat — 8n +1) (12) 


The significance of row and column effects is tested by comparing the corre- 
sponding mean square with the error mean square, SS,/(n? — 3n + 1). 


TABLE 3 
Abrasion Data! 


37 

5.03 5.50 
62 

4.96 
75 

6.34 
87 

6.31 5.46 
50 

7.27 6.54 
25 

5.93 8.02 


Analysis of Variance 


DF. SS. MS. 
Total 36 1462.8916 — 
General mean 1 1431. 1089 oa 
Rows 5 2.1897 0.438 
Column 5 2.5742 0.515 
Treatments 5 23.5300 4.706 
Error 20 3.4888 0.174 


_ Number in upper left of cell represents % relative humidity; number in low right of cell 
ig abrasion measurement. 
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APPLICATION TO A TECHNOLOGICAL PROBLEM 


The experiment mentioned in the introduction of this paper appears as an 
illustration in Bennett and Franklin’s book [3]. Table 3 contains the data and 
the analysis of variance as given by these authors. The treatments are the six 
values of relative humidity, in percent, and the measurement is rate of abrasion, 
in millimeters per 1000 machine strokes. The analysis is based on the assump- 
tion that the pattern of variability of the hide in a direction parallel to the 
backbone is the same, except for addition or subtraction of a constant, in all 
rows. By the method outlined in this paper it is possible to gain some insight 
into the validity of this assumption. To this effect, values A; are associated 
with the rows, and values B; are associated with the columns. For both A; and 
B, , we choose the set of values — 5, — 3, — 1, + 1, + 3, + 5. This set satisfies 
condition (2) and represents a linear effect in both directions. By introducing 
the term 8A,B; , allowance is made for the existence of row-column interaction, 
at least to the extent that it conforms to such a multiplicative term. The analysis, 
in accordance with Table 2, is given in Table 4. 


TABLE 4 
Analysis of Variance of Abrasion Data, Allowing for 
Row-Column Interaction 


DF. 8S. S. 
Total 36 1462 .8916 
General Mean 1 1431. 1089 
Rows 2.1897 


5 
Columns G . 2.5742 
5 


Treatments, ignoring interaction 
Row-Column interaction, corrected 
for treatments 0.6986 


23.5300 


Error 2.7902 


Row-column interaction, ignoring : 
treatments 0.9649 
Treatments, corrected for 
row-column interaction 23.2637 


It is seen that the interaction term is indeed significant at the 5% level, thus 
providing evidence that the ordinary latin square analysis may not have been 
entirely adequate. 

Bennett and Franklin pursue the analysis a little further by fitting a straight 
line to the treatment values plotted against percent relative humidity. They 
show that the mean square corresponding to the scatter of the treatment averages 
about the fitted regression line is somewhat larger than the error mean-square, 
but nevertheless fails to reach significance, even at the 10% level. The fitting 
of a straight line to the plot of treatments against relative humidity can readily 
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be carried out in the presence of a significant row-column interaction of the type 
considered in this paper. 


LINEAR Fit 1In TERMS OF A CONTROLLED VARIABLE 


To consider the problem in general, assume that instead of model (1) we 
actually have: 


Yr =e + pi +7; + BAB; + Ku — ®) + Gn (13) 


where «4, represents the controlled variable (in the example discussed here, 
relative humidity), @ is the average of the n values of u, and K is the unknown 


slope of the linear relation between treatment effects and the controlled variable. 
Let 


ti =~ i= VE (14) 
Then the least-squares estimate for K is given by: 


Oe x Piso Tw) — (Lo x YsiePisd De Pw,) 
nD b Pm) = (Lo Pw.) 


The standard error of K is: 


oR = 


Eee] 
m2) De pid) — (Lo Pan)’ | 


TABLE 5 
Analysis of Variance for Test of Linearity of Treatments Versus Controlled Variable 


Source of § Degrees of 
Variation Freedom Sum of Squares 


Row-Column [> x YiixDisl 
Interaction, SS, Li 


ignoring 7 a a Dii 
+ 7 


treatments 


I, thus ff linear ét, (XL LEME To) — (LD LX yeaa D Pwd? 
e been corrected for SS.) = ewe eee St Ss OS 

(L D pidln( 2 D pid Lo) - (Lo Pw,)'] 
traight 


interaction 
. They ff Error nt— 2-1 SS, = SS, — [SS, + SS, + SS, + SS + SS,o] 
verages Error, hee 
square, Table 1 n? — 3n + 1 S S; 


fitting Remainder n—-2 SS, = SS — SS, 
readily 
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When model (13) applies, the proper estimate for 8 is not that given by 
equation (6), but rather: 


nD » YsiePa( 2) — (Qo Tws)( 2s Pits) 
. > See) =a. — 


The standard error of this new estimate of is: 


aj Leahigaaamiacittedinesibdibeniiias’ 
yd ae Pio) — CS Par a ~ 


The relevant part of the analysis of variance for model (13) is given in Table 5. 

In order to test the adequacy of the linear fit, the mean square for error in 
Table 5, i.e., SS, is partitioned into two parts, SS; (see Table 2) and the 
remainder, SS,. . The test is carried out by considering the F ratio: 


SSi2/(n — 2) 
~ SS,/(? — 3n + 1) 


If this F ratio is found to be significant, the linear fit of trea’ 1ents versus the 
controlled variable cannot be considered to account fully ‘or the variability 
between treatments. 

Applying the foregoing to the leather abrasion data, we obtain the results 
in Table 6. 


TABLE 6 
Analysis of Variance of Abrasion Data, Fitting a Linear Trend versus Relative Humidity 


DF. 8.8. M.S. 
Row-Column Interaction, 
ignoring treatments 1 0.9649 0.965 
Linear fit vs. relative humidity, 
corrected for interaction i 21.9441 21.944 
Error 23 4.1098 .« OF 
Error, from table 3 19 2.7902 0.147 
Remainder 4 1.3196 0.330 


The remainder mean-square, tested against the error mean square obtained 
in Table 4, yields an F-value of 2.24, with 4 and 19 degrees of freedom. This 
value is less than the critical F-value at the 10% level of significance, 2.27. 
Therefore, the expression of the relation between abrasion resistance and 
relative humidity by a linear equation is consistent with the data. Thus even 
though the interaction between rows and columns is significant, it does not in 
this case change the conclusions reached under the assumption of additivity 
of rows and columns. If the linear relation of abrasion resistance versus relative 
humidity is assumed correct, the best estimate for the slope of this relation 
is given by equation (12) which, for the present data, gives a value of 

= — 0.0367. 
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ComPparRIsoN Witu OTHER Trsts For Non-AppITIvVITY 


In 1949, Tukey [4] proposed a test for non-additivity in a two-way classifi- 
cation analysis of variance without replication. Underlying this test is the 
following model, in which the notation is essentially that of equation (1): 


Yi =Beeaty+ e (po; + vi)” + €; (19) 


In a subsequent note [5], Tukey illustrated the application of his test to a 
latin-square design. 

Ward and Dick [6] have developed an iterative procedure for the analysis 
of a model of the following type: 


Yi = e+ ps +7: + Bo: + 4G: (20) 


If the row and column effects, p; and y; , represent quantitative variables, 
they may be thought of as polynomials of simpler variables, such as the A; and 
B, used in the present paper. From this point of view, the term f’p,y; represents 
a polynomial in two variables, the leading term of which is precisely the bilinear 
term BA,B; . Thus, the proposed procedure is related to the tests of Tukey 
and Ward and Dick, but restricts the interaction of rows and columns to a 
simpler term. 
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A Graphical Estimation of Mixed Weibull 
Parameters in Life-Testing of Electron Tubes 


Joun H. K. Kao 


Cornell University 


It is widely recognized that electron tube failures may be classified into two types: 
sudden and delayed. A mixture of two Weibull distributions, each representing one 
type of tube failure, is proposed, and a simple graphical method for estimating the 
parameters of the mixed Weibull distribution described. 


INTRODUCTION AND SUMMARY 


The use of Weibull distribution (1) in characterizing the lifetime of electron 
tubes was first studied by the author some years ago. Some of the results were 
given in (2, 3, 4, 5). The present paper includes a refinement in the mathematical 
model as follows. Instead of a single distribution, a mixture of two Weibull 
distributions each representing a type of tube failure is postulated. The concept 
of two failure types, sudden and delayed, concurs with the field experience 
on tube failure. This new model proposed here, allows even more flexibility in 
data fitting with very little additional complexity over the single-population 
model used earlier. A simple graphical method of estimating the parameters of 
the mixed Weibull distribution on a special graph paper is illustrated by using 
actual life-testing data gathered by a small-scaled life-test of some 800 6AQ5A’s 
conducted at Cornell University under a Signal Corps contract. The construc- 
tion of a statistical confidence band on the theoretical failure distribution 
utilizing the incomplete beta function table is also described. 


CLASSIFICATION OF FAILURE TYPES 
In broad terms electron tube failures can be classified into two types: 
Catastrophic or Sudden failures and 
Wear-out or Delayed failures. 


A tube which is operating normally and then suddenly becomes completely 
inoperative is defined here as a catastrophic failure. A tube whose character- 
istics gradually fall outside of some specifications is defined as a wear-out failure. 


1. Catastrophic or Sudden Failures. 
In general, the most important causes for catastrophic failure are of a 


Work done under U. 8. Signal Corps Contract DA 36-039-SC-64646 with Cornell Uni- 
versity. 
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mechanical kind. The catastrophic failures include glass failure, short circuits 
and open circuits. The glass failures are mainly due to strain cracks in the glass, 
and seal leaks; the shorts are due to loose conducting particles, such as lint, 
and to actual contact of parts that were spaced too close, such as pushed grid- 
wire turns; and the opens are due to open welds or burned out. elements. Often, 
there are also indirect causes. For example, the current in a tube which is a slow 
leaker may rise excessively and cause a tube element (e.g. screen grid) to melt 
and short to another element. 

For most applications, since it is usually the poor-quality tubes that fail 
prematurely, the hazard rate (to be defined later) of the population decreases 
with operation on life. Furthermore, it is reasonable to assume that failures of 
these poor-quality tubes would actually start as soon as they are exposed to 
risks. Under these considerations a Weibull distribution with location parameter 
equal to'zero and shape parameter less than one would be an appropriate math- 
ematical model. 


2. Wear-out or Delayed Failures. Ls 


The most important causes for wear-out failure are of an electrical or chemical 
nature. The wear-out failure is associated with the “‘wearing-out”’ of one or more 
of tube’s elements. In the case of a conventional receiving tube, to cite a specific 
example, the cathode element “wears-out’’ as a result of the evaporation of 
barium from the coating, the loss of coating and barium due to ion bombard- 
ment, the formation of interface resistance, and other related factors. The 
“wearing-out” of the cathode effects such important tube characteristics as 
transconductance, peak emission, plate current, ete. which often have specific 
life-test end-points. In actual field use of tubes, the slumping of tube character- 
istics causes. ‘‘wear-out’’ failure according to applications which include circuit 
requirements such as various duty cycles; environmental stress such as elevated 
temperature, maintenance policy such as marginal checking, legal regulations 
such as those imposed on airborn applications, etc. 

For most applications, very few wear-out failures occur until some particular 
time period has elaspsed after which the slumping of characteristics causes the 
hazard rate to increase with time. The length of this delay period depends upon 
application and in the case of static life-testing upon the specific definition 
of life-test end-point. In view of these considerations a Weibull distribution with 
some positive location parameter and shape parameter larger than unity seems 
to be a reasonable model. 

Combining the above two failure causes, a mixture of two Weibull populations 
as described in the following sections is postulated as the underlying distribution 
of the life time of electron tubes. 

A few words about why the Weibull distribution is used instead of, say, the 
Pearson Type III, log-normal, or other seems in order here. First of all, to use a 
statistical model to characterize a physical phenomenon such as the lifetimes 
of electron tubes whose failure mechanism is not entirely known, is at best an 
approximation. There is little theoretical justification for choosing the Weibull 


roe 
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distribution as the failure model except that its c.d.f. form coincides with the 
form of the general failure distribution derived in (5), which is not true of the 
Pearson Type III or log-normal distributions. Besides a theoretical justification, 
other criteria of a “good” statistical model seem to be: (a) It should fit the 
data, especially near the left-hand tail where the fit is of most concern; (b) It 
should be (although this is not essential) easy to handle analytically. The 
Weibull distribution is very flexible, especially the mixed model, and most derived 
characteristics of life quality can be expressed in closed form. Until a better 
model is found we choose the Weibull distribution. 


DEFINITIONS. 
We give the following definitions which are relevant to this paper. 
1. Simple Weibull Distribution. 
A simple Weibull cumulative distribution function (c.d.f.) is defined as, 
(1) F(z) =1—e'@"7"/" for 2>-7;  a,B positive. 


a = scale parameter 
8 = shape parameter 
y = location parameter* 


The first derivative of F(x) with respect to x is called a simple Weibull prob- 
ability density function (p.d.f.), thus, 


(2) dF(x)/dx = f(x) = Be = 9)P* 1-10 


a 


Figures (1) and (2) show some of the simple Weibull distributions. 


2.0 
Ficure 1—Simple Weibull c.d-f. (a = 1). 


* The notations originally used by Weibull (1) are: z» , m and 2, which correspond to our 
a, 8, and + respectively. 
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Figure 2—Simple Weibull p.d.f, (a = 1). 


2. Mixed Weibull Distribution. 
A k-fold mixed Weibull c.d.f. is defined as, 


k k 
(3) F(x) = > iF (2) where 0<p; <1, > pi =] 
{=i i=1 


and 
Fz) = 1 — exp {—(2 — y,)**/ai} 


is called the 7th sub-population in c.d.f. form. The quantities p, are the pro- 
portions of sub-population mix or simply mix parameters. If we put R;(x) = 1 — 
F,(x) then 


(4) F(x) =1—- 2d pRi(2). 


A k-fold mixed Weibull p.d.f. is therefore, 


k 
(5) f(a) = 3 patie) 
with sub-population in p.d.f. form 


fla) = BE = WD exp {—(@ — y)"/au) 


Figures 3 and 4 show a 2-fold mixed Weibull c.d.f. and its corresponding p.d_f. 
each with two sub-populations with mix parameters equal to p and 1 — p = q. 
3. Composite Weibull Distribution. 

An r-component composite Weibull c.d.f. is defined as, 


(6) F(z) = F(x) for 6; <2 bj, j= 0,1, 2, +++ 1; 
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(a, =o, p=! B a B=, p=.20, 7,70; 7p= +40) 


0.8 Ly2 1.6 
Figure 3—Mixed Weibull c.d.f. 


where 
F (2) a eo 7 rPisag 
? 


is called the jth component in c.d.f. form. The quantities 6,’s are the points of 
component partition or simply partition parameters. Since the end points are 
5) = 7, and 6,,, = ©, hence for an r-component composite Weibull distribution 
there are only (r — 1) partition parameters. 

An r-component Weibull p.d.f. is 


(7) f(x) = f(x), for 6; 2S $41, j = 0,1,2, +--+ ,7; 


(@,=a,=1, 6 \72 B B,=3, p=.20, 7)=0, 758-40) 


nae) | | | | ret | | 
Bt] | ees 


7 
Figure 4—Mixed Weibull p.d.f. 
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where the jth component is 


f;(x) we Bx aa 7°" e771 BI/ as 


a; 
Figures 5 and 6 show a 2-component composite Weibull c.d.f. and its corre- 
sponding p.d.f. partitioned at 6. *e, 


(A,=0),=1, Bz=l, 6, =3, 73=7,=0; S=1) 


i bai led | ttt 4 


0.8 


0.6 


Figure caiaiies Weibull c.d.f. 


(a sO, =1, ‘Bal, B,=5, 7377170, $=!) 


0 A 0 e 8 dX rs 2 no 26 
Fiaure 6—Composite Weibull p.d.f. 


4. Moments of Lifetime Distribution. 
Given f(z) and F(z) as the p.d.f. and c.d.f. of lifetime with y < x < @, by 
a theorem shown in Appendix B the kth moment is 


(8) wa fata eth fot — F(2)} dz 
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Equation (8) simplifies the computation of moments of Weibull distributions 


considerably. If we put 
+ 4 ‘ ™ x74 
ite v z d 
(9) [ 2” ec * dz, 
then for a simple Weibull distribution, 
(10) Mean p= y+ a°I(b + 1) 
(11) Variance o° = a” {I(2b + 1) — (b+ 1)} 


_ 1(3b +1) — 3F(2b + I)T(b + 1) + 2M) + - 1) 
(12) Skewness a; = {1(2b es 1) a rb - 1)}*? 


5. Survival Function and Reliable Life. 


Survival function R(x) which is also known as reliability function is defined as, 
(3) R(2) = [| fv) dy = 1 - FG) 


The reliabie life, p, , which is sometimes called the minimum life, is defined for 
any specified r such that, 


(4) ve r= f f(y) dy 


The reliable life p, is the same as the quantile of order p if p = 1 — r. A special 
case is when r = 3, then p, becomes the median &. In the case of a simple Weibull 
distribution, we find 


(15) Weibull R(x) = exp[—(x — )°/a] and 


(16) Weibull p, = y + a” (— Inr)’ where Inz denotes the natural logarithm 
of x. 


(17) Weibull € = y + a’ (In2) ”. 
6. Failure Rate, Hazard Rate and Retired Life. 
Failure rate for period of length h, denoted by G(z, h) is defined as, 


08) Ge, = Ef )/R@)) dy = FRE Leth) — Fe)) 


The last member of the equation is obtained by the definition of F(x). Hazard 
rate, Z(x), also known as instantaneous failure rate, is defined as, 


(19) Z(x) = Lim G(a, h) = f(x)/R(@). 
h-0 
The last member of the equation is obtained by the definition of f(x). Retired 


life, ¢, , which is sometimes called replacement life, is defined for any specified 
z such that, 


z = f(S)/R(S) 
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In the case of a simple Weibull distribution, we obtain, 


(20) Weibull @(z, h) = if ep [aaah 


(21) Weibull Z(x) = A(z — 7)*"/a 


1/(B-1) 
(22) Weibull ¢, =7 + (2) 


7. Life Expectancy. 


The life expectancy, denoted by L(x), which is the additional expected life 
on an item aged 2, is defined as, 


(23) i= / oe / ” R(t) dt/R(2) 


The last member of the equality requires that Lim,.... x” f(x) = 0 which is true 
for most distributions. 


(a) For <7, Le) = (y—-2+ [ Ro dt}/R(y) = » — x 


a mee>4, a / ” R(t) dt/R(2) = [ Rw + y) dy/R(z). 


For a simple Weibull distribution, we have for the second case, 


(24) Weibull L(z) = 5 np e—a"\" 


8. Probable Life. 


The probable life, denoted by B(x), which is the total expected life on an 
item aged 2, is defined as, 


(25) B(x) = L(x) + «. 


If can be seen that for a new item, i.e. x = 0, both L(x) and B(x) become the 
same as mean life ». Hence, Equations (23), (24), and (25) are usually used as 
measures of life quality when the items are not new. 


THE Proposep Mopet-MIxep WEIBULL DISTRIBUTION. 


The foregoing discussions on catastrophic and wear-out failure types of elec- 
tron tubes suggest the following two sub-populations: For catastrophic or sudden 
failures, 


(26) F,(xz) = 1 — exp {—2"*/a,}, for «>0O and a, >0,.0< 8, <1. 


*1r(k| a) = f° y*e-v dy which is called the incomplete gamma function and is tabulated 
in (8). 
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For wear-out or delayed failures, 
(27) F.(xz) = 1 — exp {—(x — y.)**/a2}, for <> -y. and a,>0, &>1. 
A mixture of these two sub-populations is our proposed model, thus, 
ogy Fe) = P(e) + ah 2) 
= 1 —pexp {—2"*/a,} — geap {(-(¢ - ¥2)"* /ats} 
with p + q = 1. Equation (28) states that at a given application the lifetime of 
an electron tube is a mixed chance variable depending not only on the tube’s 
mechanical strength but also on the tube’s wear-out quality. The catastrophic 
scale and shape parameters a, and 8, ; the wear-out scale and shape parameters 
a, and B, ; the delay period 2 and mixed proportion p or its equivalent g = 1 — p 
are all to be estimated from the life-testing or field failure data. The following 
graphical method of estimating these parameters is based upon the fact that a 


simple Weibull cumulative distribution becomes a straight line in In versus 
In — ln coordinates. (Refer to the example given in Table I and Figure 7). 


1. Plot the sample c.d.f. of the life-testing or field failure data on the Weibull 


probability paper and fit by inspection a curve (called Weibull plot) among 
these points.** 


2. Starting from each end of the Weibull plot draw a tangent line and denote 
them by pF , and oF, which are estimates of pF (x) and gF.(x) respectively. 


3. The intersection of qF 2 With the bottom borderline gives 7, which estimates 
Y2 


4, At the intersection of gF, with the upper borderline drop a vertical line 


whose intersection with pF, as read from the Percent Failure scale gives p which 
estimates p (and hence gq = 1 — p). 

5. Knowing the estimate of p, the mixed Weibull population may now be 
separated into the two sub-populations. The estimated number of items belong 
to F,(x) and F,(x) are np and n(1 — #) or ng respectively from which the sample 
c.d.f.’s of F(z) and F.(x) may be calculated. 

6. The y-intercept and slope of the Weibull plots of the separated sub-popu- 
lations F’;(x), give the usual estimates of Ina; and 8; respectively, (7 = 1 and 2). 


An APPROXIMATION TO THE PROPOSED MopEL-ComMposITE 
WEIBULL DISTRIBUTION. 


The above graphical procedure of estimation seems quite involved. A some- 
what simpler method is available if the mixed model can be approximated by 
a composite model. This approximation is good for small p and large 7, in the 
mixed model. In this case, the Weibull plot would consist of two distinct linear 


** For small sample size, it is better to use j/(n+1) instead of j/n for the jth sample c. d. f. 
See Appendix A for details. 

** Weibull probability paper contains In versus In-In scales and other scales calibrated 
for life-testing use. It is available from Cornell University, Ithaca, New York. 
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portions which can be treated as the two components of a composite population 
as follows: 


F(x) = F(x) = 1 — exp {—2**/as}, 
(29) for 0O<2< 65, a; > 0, 0<8 <1; 
F(a) = F(z) = 1 — exp {—2"*/au}, 
for 6<2< oo, a% > 0, B,>1. 


Equation (29) is a 2-component Weibull c.d.f. with both location parameters 
equal to zero, catastrophic scale and shape parameters equal to a; , 8; and wear- 
out scale and shape parameters equal to a, , 8, .* The partition parameter 6 at 
which the two components intersect is not a real parameter, because it can be 
expressed in terms of other parameters. At time 6 the proportion of catastrophic 
failure is equal to that of wear out failure, thus, 


(30) 1 — exp {—8**/a,} = 1 — exp {—6**/a,}, hence 


_ fox 1/(Bs—Bs) _ (io a inas) 

’ (x) biker B, — Bs 

Estimation of parameters in the composite case becomes much simpler. (Refer 
to the example given in Table I and Figure 7.) 

1. Plot the sample c.d.f. (See Appendix A for small samples) of the life-testing 
or field-failure data on the Weibull probability paper and fit by inspection two 
intersecting straight lines among the points. 

2. The intersection of these two components as read from the Failure-Age 
scale gives the estimate for 6. 

3. The y-intercept and slope of these two components gives the estimates of 
Ina; and B; , (¢ = 3 and 4). 


ConFIDENCE BaAnpD For F(z). 


Just as the Weibull plot provides a point-by-point estimate of F(x), the 
following confidence band will give an interval-by-interval estimate of F(z). 
This confidence band can be found as follows. 

In Appendix A, it is shown that for a set of ordered observations denoted 
by a, < 2% < +++ < 2, which were originally independent and identically dis- 
tributed, the probability integral transformation Y; = F(X,) for alli = 1,2,---, 
n would give the p.d.f. of Y; as the following beta distributions: 


(31) ny) = REED vid — wy" for OS we S15 


where a = tandb=n+1 —i. 


The c.d.f. of Y; denoted by H(y;) is the so-called incomplete beta function and 
is given in graphical form in Table 17 of (6) as Iy; (a, b) which is of course identical 


* The subscripts 3 and 4 used here are to distinguish from the a 1 ‘end 2 2 already 
used for the mixed model. 
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with H(y,). Knowing any three of these four quantities; J, y; , a and b, the re- 
maining one can be easily found by utilizing the following identity which was 
implied by Forms 1 and 2 of H(y,) derived in Appendix A. 


(32) I,-1,(b, @) = 1 — Iy,(a, 6) 


Tasiz I 
A Summary of Cornell Experiment, 400 Type 6AQ5A, 
Wear-Out Failure Criterion: 73.2% gm Drop. 


Estimated failure Sample c.d.f. 
lifetime in hrs. in %. 


a. Mixed population: (n = 400) 
(19.5) 
(84.4) 
(362) 
(488) 
(1050) 
1390 
1650 


ASSASSAASRSaSsE 


» say 10) 


b. Sub-population I, F(x): (np = 400 (.023) = 
1 


19.5 

84.4 
362 
488 
1050 
2390 


c. Sub-population II, F(x): (ng = 400 — np = 390, 72 = 800) 
zt t= M4 
(1390) 590 256 
(1650) 850 1.538 
(2060) 1260 3.590 
(2390) 1590 7.179 
(2650) 1850 17.949 
(3130) 2330 32.051 
(3810) 3010 55.897 
(4530) 3730 75.897 
(5060) 4260 92.564 


SVELUSE SASS oo y~-~ 


ocoanoarwhd = 


NOTES: 
1. The six points parenthesized in (a) are lifetimes of catastrophic failures, which are 
identifiable through “‘Autopsy”’ after each tube’s “death’’. 
2. For sample size as large as 400 used here, the sample c.d.f., the sample mean of F(X;) 
and the sample median of F(X;) yield almost the same value, hence only sample c.d.f.’s 
are shown. For F(z) see Table II. 
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In Appendix A, H(y;) = Iy; (t,n + 1 — 7) is equated to one half for finding the 
medians of Y,; . By the same token we can set H(y;) equal to any number between 
zero and one and obtain the so-called quantile of Y; . Hence for the confidence 


coefficient of 100 (1 — a)%, the lower confidence limit L(y;) and the upper 
confidence limit U(y;) are given by 


ne -s 

U(y) = Hy; — a@/2) for i = 1,2, --+ ,n. 
Plotting (33) with each x; we obtained the 100 (1 — a)% confidence band for 
F(a;). 

To illustrate, take n = 10, 7 = 2, the lower confidence limit of the 2-sided 
80% confidence interval for F(x.) is L(y2) = H;* (.10). This is equivalent to 
Ty.(2, 9) = .10. But this expression can not be found in the incomplete beta 
function table, hence we make use of the identity given by Equation (32) by 
stating J,_y, (9, 2) = .90 and solve for (1 — y,). This is done by referring to 


Table 17 of (6) which gives (1 —y.) = .945 or yz = .055. (See the second entry 
on the lower confidence limit in Table IT). 


Taste II 
Confidence Band (80%) for Catastrophic Sub-population, F(x) 


Time readings Time readings Est’d failure Sample Sample Lower Upper 
last inspection thisinspection _ lifetime c.d.f. Median Confidence Confidence 
in hrs. in hrs. in hrs. in % in% Limitin% Limitin % 


19.5 10 . 1.0 20 
20 . 5.5 33 
30 : 44 
40 : 55 
50 ; 64 
60 ; 73 


In Table I, the section (a) shows the abridged version of the life-testing data 
obtained from a Cornell experiment (7) on 800 Type 6AQ5A conducted under 
Signal Corps Contract DA36-039-sc—64646. 


Figure 7 shows the Weibull plot. In this case, the estimate of p = 2.3% (small) 
and y. = 800 hrs. (large), hence the catastrophic:and wear-out components 
of the failure distribution are quite pronounced. The sections (b) and (c) of 
Table I give the Weibull plots of F(z) and F.(x) which are shown in Figure 7. 
The following are estimates of the remaining parameters: 


ay Qe a3 Os, B, Be Bs Bs 6 
448 11>X10° 200 54x10° 48 33 236 4.2 14.50 


A two-sided 80% confidence band of F,(x) is given in Table II and shown in 
Figure 7 together with the median of F(z). 
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Figure 8 shows a plot of the hazard rate for the Cornell data using the estimates 
of parameters based upon the composite Weibull distribution. That is, 
362 "** 
A(z) = “900 


4.22°°? 
= 5,400,000 


for x < 1450 hrs. 
(34) 


for xz > 1450hrs. 


Some REMARKS. 


If the Cornell data are representative, then it can be seen that before the 
partition time (1450 hrs. for Cornell data), the hazard rate decreases with time 
and after that the hazard rate sharply increases. This empirical fact of first 
decreasing and then increasing hazard rate was observed by many experimenters 
on tube failures. 

The discontinuity near the partition point should be smoothed out, since the 
composite Weibull distribution is an approximation to the proposed mixed 
model. The inset in Figure 8 showed the enlarged portion of the lower left corner 
of Figure 8. Note that for a brief time-period (from approximately 500 hrs. to 
1000 hrs.) the hazard rate stabilizes to a constant value, hence for this time- 
period the failures can be considered as occurring randomly and the exponential 
distribution applies. However, this period would be shortened if the circuit 
requirement imposed on the tubes becomes more stringent. 

It is interesting to note that if the Cornell test were stopped at 500 hrs, 
as practiced by tube manufacturers on many military and commercial tubes, 
then the tubes removed from the life-testing racks would be better than the 


a 
g 
be 
8. 
we 
8 
a 
3 
i 
8 
z 


10 5 
Lifetime in 100 hrs. 


Figure 8—Hazard rate for Cornell data. 
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original new tubes. This is curious especially when such 500-hr. life-tests are 
considered as destructive tests. Hence, the widely practiced 48 hour “burn-in” 
as part of the tube manufacturing process in weeding out the so-called “weak 
sisters” seems to be in the right direction. For tubes on missile application a 
somwehat longer “burn-in” seems to be even better provided the circuit require- 
ment is not so critical as to move forward the wear-out component too much. 
For other areas, such as the industrial or computer application where tube 
longevity is important, a prolonged “burn-in” only has the effect of shortening 
the “de-bug” period and in no way improves the replacement or retired life of 
electron tubes. In fact, it consumes a portion of the otherwise useful life of the 
tube. 


AppENDIxX A 


THREE EsTIMATES OF THE F(X) PopuLation CUMULATIVE 
DIsTRIBUTION FuNcTION (c. d. f.) 


1. Sample or Empiric c.d.f. 


I. Ungrouped data: From a life-test of sample size n, let 7, < 2 < +--+ < 2, 
be the r exact failure-ages observed. The test is stopped when the rth failure 
appears where r < n. Let F(x;) be the estimate of F(z;), then the sample c.d.f., 
ie., the proportion failure on or before x; is given by: 


F(x) = ifn 


If r = n, ie., the complete data, then P(x,) = 1. Since in the Weibull plot 
the quantity Inin 1/[1 — F(z,)] is plotted as the principle ordinate, we would 
have to sacrifice the last point because In In 1/[1 — P(zx,)] = © which obviously 
can not be plotted. 

2. Grouped data: Let the failure-data be grouped into: z, , f; ; 22, fe-** ; 
2, , f, where z;, is the class-mark (a middle point) of the jth grouping cell and f; 
is the corresponding frequency of failures in the jth cell with >>*_, f; < n. Then 
the estimated c.d.f. for 7th cell is; 


F(e,) = > fi/n. 


If we have the complete data, then >—%_, f; = n and F(z,) = 1. Again, the 
last point is omitted in the Weibull plot. 

Accordingly, for complete life-testing data if the sample size is small (if the 
number of cells is small in the case where grouped data are used) the information 
represented by the largest cell is sacrificed in the Weibull plot. This disadvantage 
can be rectified if the next method of getting F(z) is used. 


II. Sample Mean of F(X;). 
Let 1, < %2 < +++ < 2, bea set of n independent observations ordered accord- 


ing to their size from a p.d.f. f(x) whose c.d.f. is F(x). Then by the multinomial 
theorem for p, + pe + ps = 1 the probability element g(z,)dz;, is: 


n! i-1, 1 ni 


G-Dlie-oim pops = C,{F(a,)}*"{1 — F(a,)}""“f(a,) de; 
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where C,, = n!/(¢ — 1)! (n — 1)!. Allowing a change of variable y, = F(z;) or 
dy; = f(x,;)dx; the p.d.f. of y, is, 


hty;) = Ciyi (1 = yi)" for O< y, <1, 


since F(z,) = F{F"‘(y,)} = y; and h(y;) = g{F7*(y,)}. |J| where J is the 
Jacobian of the transformation and, 


J = dz;/dy; = 1/{dF(z;)/dx;} = 1/f(x,). 


Since n and 7 are integers, we may write C, as T'(n + 1)/{T(@) T(n + 1 — 2).} 
The density h(y) can now be recognized as a beta distribution with parameters: 
a = iand b = n + 1 — i, where the beta distribution (p.d.f.) is defined as: 


fea) = Te ara - ay fr OSekt a b30 


From the identity that a! b! = (a + b + 1)! fj 2°(1 — x)*dz, the kth moment 
of y; about the origin is, 


n! (K+7—)!(n—7! 


Ey) = C, [ gL — yy" dy = i Hie ar (+n)! 


_ al(e+i—D! 
~~ G-DIk+n)! 
Hence 


ni 1! a 


EW) = G@o)l@+ Dl atl 


Recall that y, = F(x;), the expression above states that 7/(n + 1) is an unbiased 
estimate of F(x;) hence for, 


1. Ungrouped data: The sample mean of F(X;,) is given as: 
F(a) = i/(n +1) 


If i = n, i.e., the complete data, then F(a.) = n/(n + 1) which is still within 
graph for n less than 1000. 

2. Grouped data: In this case, the sample mean of F(X,) is given by: 
F(z:) = di-. f;/(a + 1), where 2; for each 7, are defined as before. Again, if 
we have the complete data, then >-*., f; = n and F(z) = n/(n +1) = 1-—- 
1/(n + 1) and In In i/{1 — F(z,)] = In In (n + 1) which is finite for any finite 
sample size n. 


III. Sample Median of F(X;): 


From h(y;), the p.d.f. of y; derived above the c.d.f. of y; , denoted by H(y;), 
may be found. The median of y; is simply the solution of H(y;) = 3. The c.d.f. 
of y; is, 
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Hg) = C, [20-2 de = 0, fF (0! (05) ae 
= 6. eo (G 1) [eran oe Zen G yao 


(Form 1) 
If we let v = 1 — z, the above integration may be evaluated as: 


Hiy) = C, f (1 — »)'“v"-! dy = C, . oF (1)! (' 7 " v! dv 


i=0 


=¢, 5 (-n! (' 5 " eo dy 


i=0 


‘ia ¢. > (-0 ( j ') {1 ex (1 a yg +n fe ae i)? 


7=0 


(Form 2) 
Using Form 2 of H(y), the setting of i = 1 gives, 


Hy.) = 1— (lL — m)", 


H(y:) = 4 gives y, = 1 — (3)'" = F(a) 
Using Form 1 of H(y), the setting of 7 = n gives, 


H(y,) = yh 
hence H(y,) = } gives y = (3)'” = F(a,) 


These two explicit solutions give the estimates of F(X;) by its median at 2, 
and x, . For other values 22 , 23 + X,-; (as a matter of fact x, and x, also), the 
median may be found by setting H(y,;) = 3 in Table 17 of (6) and solve for y; . 


AppENpDIx B 
Moments or A Mixep WEIBULL DISTRIBUTION 


In order to find the population mean, variance, etc. of a mixed Weibull dis- 
tribution, we shall prove a theorem which simplifies moment computations for 
many distributions for which 1 — F(z) is in the closed form. 

Theorem: Let X be a random variable with c.d.f. = F(x) and p.d.f. = f(z) 
and reliability function, R(x) = 1 — F(z) all defined over y < x < © wherey 
is any real number, then, 


al 3 [ hide at eh / * R(x) dar 


if Lim,... z*** f(z) = 0 for any k < &. (note) 
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Proof: Rewrite the last member of the above expression for yu; as follows: 


/ "11 — F(a) jhe! de 


Denote u = 1 — F(x) and dv = k x*~'dz and integrate this expression by parts. 
We get, 


a*{1 — + ui 
7 


The upper limit of the first term vanishes if Lim,... z**’ f(x) = 0 for any finite 
k and F(y) = 0 by definition of c.d.f. Hence, the theorem is proved. 
Corollary 1. If the distribution is defined over some finite range y < x < 6, 
then F(y) = 0 F(8) = 1, the condition Lim x*** f(x) = 0 is unnecessary. 
Corollary 2. If the distribution is defined over the range 0 < x < o, then 


=k [2 R@) de. 
0 
As an example we can use the above theorem for finding the moments of a 


mixed Weibull distribution where Lim,... x**’ f(z) = 0. The kth moment about 
the origin is, 


w= pk [ x” exp {—2"*/a,} dx + qvi + ak [ xz* exp {—(x — ¥2)"*/a2} dx 
gi b 


We proceed to evaluate the last term by letting y = (x — 72)**/a, . Dropping 
the subscript “2” of parameters for ease of typing, the last term of yj is, 


ok [2 —(z-7)8/a dz - af (ai /Fy 1/8 + yeep ty" -1 dy 


= gk : (' - 1) Agito" er(é mae 


The first term in z{ is found similarly. Hence, 
k-1 = : , ae 
pi se pay" T r(4 + 1) fe ay: fe gk 2 r j Ntay-07 e'r(é ma ‘) 


By setting k = 1 in the above expression, we get 


= n= pol’ 1(L + 1) + dr ~ air (E ae 1) 


which is in fact the weighted average of the means of the sub-populations in 
accordance to the mixed parameters p and q. The higher moments, such as 


Note: For y = 0, k = 1, we have » = {® R(x)dz which is well known in actuarial statistics 
for the mean of a non-negative random variable (see, for example, Kurtz, E. B. Life Expectancy 
of Physical Properties). In actuarial case, since R(x) is aaaee in range in both directions the 
condition Lim—22f(x) = 0 is automatically satisfied. 
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variance and skewness, have more complicated relationships although tedious 
but can be evaluated similarly. 
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Evaluation of Chemical Analyses on Two Rocks 


W. J. YoupEN 
National Bureau of Standards 


Interlaboratory and round robin test programs can provide information, not 
only on a test procedure under investigation, but on the participating laboratories 
as well. A simple graphical technique is proposed to aid in the comparisons between 
laboratories. 


INTRODUCTION 


Test methods in industry are often submitted to an inter-laboratory or round 
robin test program. The avowed intention is to find out how satisfactory the 
test procedure is in actual practice. The point of view may be taken that the 
participating laboratories are just as much under study as the test method. If 
eight or nine laboratories out of ten get good results, perhaps attention should 
be directed to the laboratories having difficulty instead of to the test method. 
A published example of an extensive and carefully controlled study of the 
chemical analysis of rocks is reviewed in this paper with particular consideration 
given to the participating laboratories. 

Geological Survey Bulletin 980, ‘A Cooperative Investigation of Precision 
and Accuracy in Chemical, Spectrochemical and Modal Analysis of Silicate 
Rocks” 1951, reports chemical analyses by 34 cooperating laboratories on a 
granite rock and a dibase rock. The analytical results from the various laboratories 
posed a difficult problem of evaluation. The report discusses the use of the modal 
value, the use of the mean of all of the results, and the use of a mean based upon 
a consensus of more or less closely clustered values. The report includes scatter 
diagrams of the results for each determination. A separate scatter diagram was 
made for each rock. The selection of the results to be included in any consensus 
cluster is not an easy one. This note approaches the problem by an extension of 
the familiar scatter diagram. The purpose is to pick a subset of the laboratories 
and use the results from these laboratories for evaluating the analytical results. 


ANALYsIS OF DaTa 


Thirty of the laboratories reported on both rocks. These data are recorded 
in Tables 1 and 2 of the Geological Survey Bulletin 980. Each of these thirty 
laboratories makes available a pair of results; one for each rock. Consider the 
SiO, determinations. Plot these paired values using the z axis for the granite 
results and the y axis for the dibase results. Each laboratory provides a point in 
a two dimensional scatter diagram shown in Figure 1. In theory, for statistically 
independent measurements with no systematic differences between laboratories, 


409 





J. YOUDEN 


re) ° ° 
= . > 5 
= os - ° = 


Ficure 1. Dibase rock results (y-axis) plotted against results obtained on granite rock (x-axis). 
The numbers in the graphs identify the outlying laboratories. 
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TaB1E 1. Tabulation of laboratories with very extreme resulis* for siz of the elements. 


Laboratory Number 


Element 


ma ee Me A Ut 8 23 31 33 


« The asterisk denotes very extreme results. The 16 laboratories not listed had no extreme 
results for these six elements. 


the points should be scattered in a circular pattern centered on the center of 
gravity of the whole group. Almost always the pattern of points departs from 
theory in two noticeable ways. First, there is a main body of points, more or 
less elliptical in shape with the long axis of the ellipse inclined at an angle of 
approximately 45° to the z-axis. Second, a small minority of the points are 
individually separate from the main cluster and, in most cases, these points 
tend to lie near the extended long axis of the ellipse. 

Both these departures from theory have a reasonable explanation. The center 
of gravity of the main cluster represents the consensus of the laboratories for 
the method employed in the determination. The coordinates of this centroid 
are probably the best estimates of the unknown values that can be made from 
the data. The displacement of this centroid, from the true values for the two 
rocks, is largely along a 45° line through the unknown true value because any 
systematic error in the procedure may be considered to apply (in most cases) 
equally to both rocks. If the method, as used, tends to give high results by P 
percent, it presumably does this for both materials and thus the x and y dis- 
placements are equal and the point is displaced along a 45° line. Now any one 
laboratory is apt to have its own systematic error relative to the systematic 
error in the procedure itself. For the moment forget about the true values and 
focus on the centroid which represents a working datum of reference. A labora- 
tory which tends to get results higher than the other laboratories will be dis- 
placed from the centroid upward along this 45° line because it gives the same 
positive bias to each material. A laboratory with an excessively large relative 
bias will be considerably displaced along the 45° line. The general elliptical 
shape results from the collection of these individual biases. 

These two dimensional scatter diagrams are useful in the evaluation of the 
results. There is a reinforcement of the judgment in setting aside certain de- 
terminations when the scatter diagram shows an emphatic discrepancy by one 
laboratory for both rocks. There are a few instances where the displacement 
from the main group is accounted for by one of the pair of results and one cannot 
help but wonder if some simple slip in calculating or typing is the explanation. 
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In any event an examination of the scatter diagrams for SiO, , CaO, K.0, Na.0, 
MgO and Al.O; , shown in Figure 1, shows approximately six points in each 
diagram well removed from the main group. The two dimensional departure 
of the points from the centroid gives a good basis for suspecting these out of 
line results. 

The foregoing does not exhaust the inferences to be drawn from these scatter 
diagrams. The elliptical character of the cluster has been pointed out. The 
more vulnerable the analytical procedure is to individual laboratory bias, the 
greater the tendency of the main cluster to assume an elliptical pattern. By this 
method of plotting a visual comparative appraisal of different procedures (for 
the same or different elements) may be made. Procedures that give a generally 
circular cluster indicate that, should a systematic error be present in the pro- 
cedure, all the laboratories are afflicted with this same systematic error both in 
direction and magnitude or that the systematic errors are small with respect 
to the precision of the method. 


TaB Le 2. Tabulation of Statistics Calculated From the Data. 


Table 20, Bull. 980 Accepted laboratories 

Element - 

and Av. of Av. of Av. 8.D2 §8.D.., 

rock all consensus 
Si0. 72.22 72.45 
52.25 52.45 
Al.0; 14.44 14.30 
15.23 15.10 
Mg0 0.39 .0.45 
6.52 6.65 
Cad 1.42 1.35 
10.95 10.95 


Na0 3.26 
2.05 


K:0 


G 
D 
G 
D 
G 
D 
G 
D 
G 
D 
G 
D 


5.51 ‘ j 
0.71 .063 .014 


® Rejection of all the data from laboratories 16, 18, 4, 15, 22, 26, 11, 14 and 21. Also the 
determinations with an asterisk for laboratories 1, 8, 23, 31, and 33 in Table 1. 

b The precision component, i.e., the standard deviation of a single determination. The 
next column gives the standard deviation for the average of the accepted laboratories. These 
standard deviations are based on an estimate of the precision. 

¢ F is the ratio of the square of the standard deviations for laboratories and precision. A 
ratio close to unity is expected if the laboratories all have the same systematic error. Values 
of F larger than about two are considered good evidence that the laboratories have individual 
systematic errors that differ from laboratory to laboratory. An F value of 3 indicates a syste- 
matic error component about equal to the precision component. 
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The collection of scatter diagrams makes possible a reinforcement of the 
rejection procedure. The laboratories identified with the outlying points for the 
determinations of SiO, , Al,O; , MgO, CaO, Na.O, and K,0 are shown in Table 1. 
These points are also identified in Figure 1. All graphs use the same unit apace 
for one percent to facilitate comparison of the procedures. 

Certain laboratories are very much in evidence as the source of the outlying 
points in these six diagrams. In fact six of the 30 laboratories account for 25 
of the 39 outlying points. The average number of outlying points for all 6 tests 
combined is about 1.3 per laboratory. If the outlying location is only a matter 
of chance, we may use the Poisson distribution to calculate the expected pro- 
portion of laboratories with none, one, two, or more outlying points. We could 
expect just about one of the 30 laboratories to be out in four or more of the six 
diagrams. Instead, six laboratories achieve this unenviable distinction. Two of 
these laboratories are remote in five of the diagrams. We would also expect 
only eight or nine laboratories never to be out of the pattern if this event de- 
pended on chance. It is encouraging to find sixteen of the thirty laboratories 
avoid the outfield. Five of the laboratories turn up in the outfield just once, 
and in two of these cases just one of the rock determinations is responsible. 
Except for the particular determinations in question, the work of these five 
laboratories appears appropriate for retention. This leaves 21 or 20 laboratories 
for determining a consensus both as to mean and dispersion for the analytical 
procedures. 

Scatter diagrams for Fe.O; and FeO were also prepared. The one for FeO; 
showed that considerable difficulty was experienced with this determination 
because about half of the points showed extremely large departures from the 
residual compact group. Almost all the laboratories that had trouble with other 
determinations obtained poor results with the Fe,O0, determination. Similarly 
the eight outlying points for the FeO results all came from laboratories having 
trouble with three or more of the six determinations shown in Figure 1. The 
FeO data present a special problem because the dibase results show a much 
greater dispersion than the results for the granite. 


RANDOM AND SYSTEMATIC ERRORS 


The two dimensional scatter diagrams illustrate an important partition of 
the errors of the reported results. There are first the purely random errors which 
are revealed by the scatter of repeated determinations made on the same rock 
within a given laboratory and second the individual laboratory systematic 
errors measured from the consensus of the results from the qualifying laboratories. 
Earlier it was pointed out that if there were no differences in the systematic 
errors, the pattern of points would tend to be circular. On the other hand, if 
the individual laboratories had perfect precision, i.e., in this case could check 
themselves exactly to 0.01 percent, the scatter would arise solely from the indi- 
vidual systematic errors. Indeed, in this improbable state of affairs, the points 
would lie exactly on the 45° line through the centroid. Inasmuch as both types 
of error are present, the points depart more or less from the 45° line but retain 
this line as the major axis for the scatter of the points. 
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Fiaure 2. Three ways of plotting results for ferrous oxide on dibase rock (y-axis) and granite 
rock (z-axis). The numbers in the graphs identify the outlying laboratories. 
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The departures of the points from the 45° line, as measured by the lengths 
of the perpendiculars from the points to the 45° line, correspond to the precision 
component of error. The lengths of the segments along the 45° line from the 
centroid to the feet of these perpendiculars measure the relative systematic 
errors. It is instructive to square each length along the line and sum these squares. 
Also sum the squares of all the perpendicular lengths. The combined total of 
these two sums will be found equal to the sum of the squares of the deviations 
from the mean for the granite added to a similar sum for the dibase. The sum of 
the squares of the deviations from the mean is a familiar step in calculating a 
measure of the dispersion of a set of results. Given just one set of results, say 
those for granite, no separation between precision and systematic components 
of errors is possible. When the results on two materials are available, the position 
of a point representing the two results may be represented by its distance (from 
the centroid) along the 45 degree line and by its perpendicular distance from 
the line. These two distances are identified with the relative systematic error 
and the precision error respectively. Given results on two materials, the disper- 
sion estimates can be expressed in this alternative and more revealing form. 
Some readers may wish to consult a statistical text on the analysis of variance 
of data with a two way classification. 

It must be pointed out that when rocks of sharply different composition are 
involved (either as to amounts of the elements or the presence of possibly 
troublesome elements) the precision component will also include any change 
in the systematic errors for the laboratory. If the two materials are rather similar, 
there is every reason to expect the systematic error to be the same for both 
materials in any one laboratory. When the range of the scatter of the results is 
about the same for both rocks, this complication is not likely to be present to 
any important degree. Another complication is the possibility that the precision 
error depends upon the amount of element present. Thus if there is a ten-fold 
difference in amount of element present in the two rocks, the scatter along one 
axis may be visibly greater along one axis than along the other axis and this 
in turn will cause a departure of the major axis of the ellipse from the expected 
45°, In this event, it is often helpful to convert the data to logarithms before 
plotting or undertaking a statistical analysis. This statistical device will give 
equality of variance under the assumption of a constant coefficient of variation. 
The FeO determinations afford an excellent illustration of this complication. 

In Bulletin 980 there is a notation that some of the results listed are based on 
duplicate determinations. No allowance was made for these duplicates because 
experience indicates that the systematic errors dominate the precision com- 
ponent. The effect is to bring these points a little closer to the 45° line than 
single determinations would fall. In consequence the estimate of the precision 
has been slightly favored. The other assumption made is that of a common 
precision for all the accepted laboratories. This is a reasonable approximation 
because the differences in precision are minor compared with differences in 
systematic error and, in any event, many determinations would be required 
to establish that such differences in precision existed. Observe also that many 
of the points far removed from the centroid are reasonably close to the 45° line 
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indicating that even the presence of large systematic errors is not necessarily 
tied in with precision. The estimate of precision, using the data from the selected 
laboratories, provides a basis for predicting the extent of agreement that can 
be expected between laboratories that in one way or another adjust their pro- 
cedures by reference to a standard sample. 

The scatter diagram first prepared for FeO showed a very much greater spread 
for the dibase results than for the granite rock. (Fig. 2). Consequently the 
logarithms of the data were used to prepare a graph. After this transformation 
there was a greater spread in the granite values than for the dibase values. 
There are two rather well defined clusters in this diagram and this suggests 
the important clue that there were two different analytical procedures. The 
centroids for the two clusters differ mainly in the means for the granite rock 
suggesting that it is on this rock that the two analytical procedures disagree. 
Finally the data were plotted using the square roots of the data. The axis of 
the major cluster is now much closer to making the expected 45° angle with the 
x axis. The suggestion of two distinct groups in the pattern continues in the 
diagram based on the square roots. Statistical estimates based on composites 
of two differing groups have no chemical utility. The examination of these data 
puts the problem of FeO determination right back in the laboratory. 

The two dimensional scatter diagram shows the need for a careful consider- 
ation of the method for determining FeO. The 45° line through the centroid of 
the data should run through the mass of points in such a way that a ruler pushed 
along perpendicular to the 45° line should encounter points on the left and 
right of the line in a random order. Any runs of points on one side, succeedéd 
by long runs on the other side are a signal for caution. Statements about the 
procedure that hold only for a particular rock are not helpful. Bulletin 980 did 
segregate the data on the basis of the procedure employed wherever this infor- 
mation was available. Not much can be learned from a result unless it is ac- 
companied by information on the method used. 

In Table 2 there is shown the precision component for the determinations of 
SiO, , Al,O; , MgO, CaO, Na.O, and K;0O. This standard deviation shows the 
optimum performance that could be expected for these determinations should 
the laboratories achieve the same systematic error. As a matter of fact, these 
very rocks, if the compositions become accepted, are intended to serve as 
standards and would, by providing a standard reference point, achieve this goal 
of a common systematic error. The amount of this common systematic error is, 
of course, exactly the amount by which the designated composition departs 
from the unknown true composition. 

Besides the precision component there are listed the averages of the results 
from the selected two thirds of all the laboratories. It is interesting to compare 
this column with both the averages for all laboratories and the consensus that 
have been taken from Table 20 in Bulletin 980. The precision component does 
not contain: the component of error arising from variation in the systematic 
error among the laboratories. A measure of this source of error is provided by 
the column of F values. Large values of F show a real need for standardization 
among the laboratories. 
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In Bulletin 980 emphasis is placed (p. 12) upon the possibility that some 
extreme value may in fact be closer to the true value than an average based 
upon all the data, a consensus of selected results, or some modal value. This 
invites the reply that whatever comfort this may be to a laboratory with a point 
at one end of the 45° line, the same gesture adds to the discomforture of those 
laboratories on the opposite end by making them even further off from the true 
value. No one will seriously contend that the centroid is without bias. But in a 
practical world, men must adhere to some standard, even a standard made by 
men. The laboratory with an extreme position along the 45 degree line may 
be the most nearly ¢orrect laboratory and the other laboratories all afflicted 
with a systematic error of the same sign, but evidence is required to support 
such a position. Unless, and until, such evidence is forthcoming, the use of a 
consensus value will persist. The double sample plotting scheme described in 
this paper should prove useful in selecting those laboratories doing consistent 
work near the consensus. 
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Errata 


“The Application of Random Balance Designs” 
T. A. Bupne, TECHNOMETRICS, Vol. 1, No. 2, 1959 


Page 141, 7th line from bottom: “+10” should read “‘—10” 

Page 141, 4th line from bottom: “B” should read “C”’ 

Page 143, Figure II: Test Run No. 5 under column headed A should be + 

Test Runs 22 and 23 under column headed G should have signs interchanged. 

Test Run No. 12 under columns headed Y; , Y, , and Y; should read 98, 93 
and 98, respectively, to be consistent with other numbers in the table. However, 
the published values of 97, 92 and 97, respectively, are the ones which have 
been used in the analyses and calculations. ; 

Test Run No. 24 under column headed Y, should read 97 instead of 102. 
This error is typographical and all ensuing calculations use the correct value 
of 97. 

The text indicates corrections of 8 units for the JK interaction and 4 units 
for the C effect. Actually, corrections of 9 and 5 units, respectively, were made 
on the data. Based on their respective analyses, the 8 and 4 are correct. Such 
a revision would make no appreciable change in the final results. 


“Response Surface Designs for Three Factors at Three Levels” 
Rosert M. DeBaun, TECHNOMETRICS, Vol. 1, No. 1, 1959 


Table 1, which appears on page 5, contains several errors. The corrected 
table is given below. 


p=1/2 p = 1.00 p = 1.50 p = 2.00 


Face Edge Corner Face Edge Corner Face Edge Corner Face Edge Corner 


1) 3° Factorial 

Design .162 .165 . -143 .188 .211 .051 .090 .121 .016. 
2) Cube + 

Octahedron .236 .248 . 1D. .253 .034 .076 .1382 .011 . 
3) Cube + 

Octahedron .272 .289 . : : . 3 .072 .124 
4) Cube + Two 

Octahedra ; ‘ ; ; ‘ P .042 .074 . 
5) Cube + Two 

Octahedra P q : ‘ ‘ 2 .038 .068 . 
6) Cuboctahedron . i ‘ 3 3 ; .062 .083 . 
7) Cuboctahedron . : P : F : .058 .077 . 
8) Cuboctahedron . ‘ : ’ ; : é .068 . 
9) Cube + 

Cuboctahedron . ‘ : f ? ‘ -044 .090 . 
10) Cube + 

Cuboctahedron . , ; ; ‘ 4 .041 .075 . 
11) Cuboctahedron 

+ Octahedron . 4 4 i i ; -056 .063 . 
12) Cuboctahedron 

+ Octahedron . ; ‘ P r ‘ .051 .058 . 





TECHNOMETRICS Novemser, 1959 


Acknowledgment 


The editors of Technometrics wish to acknowledge the editorial assistance 
of the following individuals: 


William R. Allen, H. P. Andrews, F. J. Anscombe, D. W. Behnken, E. Bianco, 
R. 8. Bingham, Allan Birnbaum, Edgar M. Bisgyer, G. E. P. Box, Irving Burr, 
Rajan Chanmugan, D. G. Chapman, W. E. Deming, Harold Dodge, A. J. Duncan, 
Benjamin Epstein, D. A. Gardiner, D. P. Gaver, Jr., R. Gnanadesikan, A. H. E. 
Grandage, B. G. Greenberg, F. E. Grubbs, Max Halperin, H. O. Hartley, M. J. R. 
Healy, David C. Hurst, G. M. Jenkins, J. H. Kao, R. W. Kennard, T. Koehler, 
Gerald J. Lieberman, H. L. Lucas, Colin Mallows, Paul Meier, Jack Moshman, 
Mervin Muller, Paul Olmstead, Emanuel Parzan, Roger Pinkham, Maynard 8S. 
Renner, F. S. Riordan, Joan Raup Rosenblatt, Judah I. Rosenblatt, Henry 
Scheffé, Willis A. Shell, Jr., C. T. Shewell, Harry Smith, Jr., Milton Sobel, 
Louis Tanner, M. E. Terry, John Tischendorf, John Tukey, H. Van Blarigan, 
G. S. Watson, Mason E. Wescott, J. C. Whitwell, S. S. Wilks, W. H. Williams, 
C. B. Winsten, William Wolman, J. Youden and William P. Youngclaus, Jr. 

The editors also wish to extend their warm thanks to all those who have con- 
tributed both directly and indirectly to the 1959 issue of T’echnometrics. 





TECHNOMETRICS Novemser, 1959 


NOTICES 


Preliminary draft of the program of the 119th Annual Meeting American 
Statistical Association of special interest to members of the Section of Physical 
and Engineering Sciences, to be held in Washington, D. C., Dec. 27-30, 1959 





Sunday, December 27—10:00 A.M. 


MULTIVARIATE ANALYsIs I 


Some Applications of Multivariate Chebyshev Inequalities: I. Olkin 
Multivariate Correlation Models with Mixed Discrete and Continuous 
Variables: R. F. Tate 


Multivariate Chebyshev Type Inequalities: A. W. Marshall 


STATISTICS IN SOME ASPECTS OF METEOROLOGY 


Problems and Progress in Weather Modification Experiments: Glenn Brier and 
Herb Thom 


Discussants: O. Kempthorne and Max Woodbury 





Sunday, December 27—4:00 P.M. 


DESIGN OF EXPERIMENTS II 


Designs of Type B’: George Box 
Paper to be announced: John Mandel 
Discussant: H. O. Hartley 


Monday, December 28—10:00 A.M. 


ScREENING Data AND Direct SEARCH BY MEANS OF 
ELECTRONIC COMPUTERS 


A Screening Program for the IBM 704: Milton Terry 
Direct Search: Robert Hooke 
Discussants: Jack Moshman and Fred Mosteller 





422 
Monday, December 28—2:00 P.M. 


MUouLttTiIvariaTe ANAtysis II 


A Representation for Anderson’s W Statistic: A. Bowker 

An Asympotic Expansion for the Distribution Function of Anderson’s Classifi- 
cation Statistic W: R. Sitgeaves 

A Representation of the Wishart Matrix: R. J. Wijsman 


Monday, December 28—2:00 P.M. 


CoNTRIBUTED PAPERS 
Estimation in Non-Linear Least Squares Problems: H. O. Hartley 
Truncation and Tests of Hypothesis, II: Irwin Guttman 
On a Significance Test of System Reliability: David Rubenstein 


Analysis of Survival Data by Regression Techniques: Scott Krane and R. J. 
Buehler 


Monday, December 28—4:00 P.M. 


MULTIVARIATE ANALYSIS, CLASSIFICATION AND DISCRIMINATION 

Classification into Multivariate Normal Distributions with ae Covariance 
Matrices: T. W. Anderson 

Some Problems in Multivariate Classification Using Discrete Variates: W. G. 
Cochran and C. E. Hopkins 

A Linear Multivariate Model for Paramorphic Repenentiiion of Clinical 
Judgment: P. J. Hoffman 

Analysis of Multi-Factor Multi-Response Experiments with Mixed Factors and 
Response Types: 8. N. Roy . 


Tuesday, December 29—10:30 A.M. 


SraTIsTicAL ANALYSIS OF TURBULENCE ON THE SURFACE or THE SUN 


Papers by: Francois Frenkiel and Martin Schwarzchild 
Discussants: G. Watson and Leo Tick 


Wednesday, December 30—8:30 A.M. 


RELIABILITY 


A Markov Model: William Allen 

On Drawing Inferences in the Case of an Exponential Distribution with 
Censoring: Jack Nadler 

Discussants: Cuthbert Daniel and Milton Sobel 
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FourtH TECHNICAL CONFERENCE OF CHEMICAL Division ASQC 


The 4th Technical Conference of Chemical Division ASQC will be held in 
Chicago on November 3 and 4, 1960. Further information concerning the 
meeting can be obtained from 


Karl J. Baue 
Abbot Laboratories 
North Chicago, Illinois 


Copies of the Transactions of the 3rd Annual Technical Conference of the 
Chemical Division, ASQC can be obtained for $2.50 from 


Jerrold H. Moyer 

Supervisor, Statistical Services 
Champion Paper & Fibre Co. 
P.O. Box 872 

Pasadena, Texss 


A list of the speakers, and papers, appears in the August issue of Technometrics. 


NATIONAL BUREAU OF STANDARDS 


PosTDOCTORAL RESEARCH ASSOCIATESHIPS 


Research associateships, supported by the National Bureau of Standards 
and awarded on recommendations of the National Academy of Sciences— 
National Research Council, are offered to provide young investigators of unusual 
ability and promise the opportunity for basic research in various branches of the 
physical and mathematical sciences. It is expected that approximately 20 
awards may be made in a total of twenty-nine fields, of which the following 
are of particular interest to statisticians: Pure and Applied Mathematics, 
Applied Mathematical Statistics, Numerical Analysis. These research asso- 
ciateships are open only to citizens of the United States, and in the foregoing 
fields are tenable only at the National Bureau of Standards in Washington, D. C. 
Applicants must have received or completed the requirements for a Ph.D. or 
Sc.D. degree, or equivalent, in one of the fields listed above at the time of entering 
upon the research associateship. 

The annual gross stipend will be $7510 and will be subject to income tax. 
Travel and moving expenses of the Research Associate and his family from 
place of residence to Washington, D. C. will be paid by the National Bureau 
of Standards. Awards will be made about April 1, 1960. Unless otherwise ar- 
ranged the tenure of a research associateship may begin after July 1, 1960 
and continue for one year, with provision for a vacation period. 

Requests for application forms or for additional information should be ad- 
dressed to the Fellowship Office, National Academy of Sciences-National Re- 
search Council, 2101 Constitution Avenue, N. W., Washington 25, D. C. Appli- 
cations for the academic year 1960-1961 must be received in the Fellowship Office 
no later than February 1, 1960. 
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Tue Rutcers UNIversity Statistics CENTER 


Rutgers University has established a Statistics Center which will be the 
State University’s central unit for research and teaching in the field of statistics. 
It will continue to offer at Rutgers the Master’s Degree program and in addition 
will have responsibility for a program of study leading to the Ph.D. in Applied 
and Mathematical Statistics. 

The Center, which is under the administration of Dr. Marion A. Johnson, 
dean of the Graduate School, will be directed by Dr. Ellis R. Ott, who has 
previously served as professor of Mathematics and Chairman of the Mathe- 
matics Department of University College, as well as Chairman of the Rutgers 
program in Applied and Mathematical Statistics. Director of Research at the 
Center will be Dr. Martin B. Wilk. The staff will also include Dr. Roger S. 
Pinkham, Dr. Mason E. Wescott and Harold F. Dodge. 

The Statistics Center will be responsible for programs of research in statistical 
theory and methodology and in ways of promoting the more effective use of 
statistics in science, industry, education and other areas. 

The Center will also make available, both within and without the University, 
consultative services with respect to the use of statistics and statistical tech- 
niques in research or other experiments and surveys, and with respect to the 
analysis and presentation of data. 

Since 1952, Rutgers has had course work leading to a Master’s degree in the 
field of Applied and Mathematical Statistics. Through June, 1959, a total of 
85 Master’s degrees in applied and mathematical statistics had been awarded, 
mostly to full-time employees in nearby industries. The present development 
constitutes a major extension of the Rutgers Statistics Program plus an ad- 
ministrative reorganization. 

The establishment of the Statistics Center is predicated on the belief that a 
group within the University charged with instructional, research, and consulta- 
tive functions in the field will both advance and strengthen the University’s 
program in statistics, and also further the effective use of statistica! tools and 
techniques by the many members of the University’s staff who, though not 
specialists in statistics, must of necessity use statistics in their own work. 
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Research Laboratory 
STATISTICIAN 


e A national consumer products company located in the midwest 
has an unusual opportunity for a statistician with administrative 
abilities to work with our research laboratory personnel. 
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.. Electronic computer available for data processing, but computer 
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PREPARATION OF MANUSCRIPTS 


Manuscripts should be submitted to the office of the editor: J. S. Hunter, 
Mathematics Research Center, U. 8. Army; The University of Wisconsin; 
Madison, Wisconsin. Each manuscript should be typewritten, double s-aced, 
with wide margins at sides, top, and bottom. The original should be submitted 
with two additional copies, on paper that will take corrections. Dittoed or 
mimeographed papers are acceptable only if completely legible. Footnotes 
should be avoided and reple.ced by remarks in the text, or placed in an appendix. 
Preferably, references in the manuscript should appear as (Jones, A. B., 1958), 
and again later in alphabetical order in a list of references. Alternatively, refer- 
ences may be numbered, e.g. [1], as they appear in the manuscript and be listed 
in this sequence in the list of references. In the reference list, each reference 
should contain, in the order indicated, the name and initials of the author 
followed by those of the co-authors, date of publication, title of reference, 
source, volume number and page. References to books should include pub- 
lisher’s name and location. 

Figures, charts, and diagrams should be professicnally drawn on plain white 
paper or tracing cloth in black India ink twice the size they are to be printed. 
A full page diagram, in print, measures 7.25 X 4.75 inches. 

As far as possible, formulas should be typewritten and symbols not available 
on a typewriter carefully inserted in ink. Authors are asked to keep in mind the 
typographical difficulties of complicated mathematical formulae. The difference 
between capital and lower-case letters should be clearly shown; care should be 
taken to avoid confusion between such pairs as zero and the letter O, the numeral 
1 and the letter J, numeral 1 used as superscript and prime (’), alpha and a, kappa 
and k, mu and wu, nu and », eta and n, etc. Subscripts or superscripts should be 
clearly below or above the line. Bars above groups of letters (e.g., log 2) and 
underlined letters (e.g., x) are difficult to print and should be avoided. Symbols 
are automatically italicized by the printer and should not be underlined on 
manuscripts. Boldface letters may be indicated by underlining with a wavy line 
on the manuscript; boldface subscripts and superscripts are not available. 
Complicated exponentials should be represented with the symbol exp particu- 
larly when appearing in the text, that is, 

exp [(a” + b”)’””] should be used in place of e*"**”"”. 
In writing square roots the fractional exponent is preferable to the radical sign. 
Fractions in the body of the text (and when possible in displayed expressions) 


and fractions occurring in the numerators or denominators of fractions are 
preferably written with the solidus; thus 


(a + D/le + d) rather than >. 


Authors will ordinarily receive only galley proofs. Fifty offprints without 
covers will be furnished free. Costs for additional reprints and covers can be 
furnished on request. 
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